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^ The phage A infection of an E. coli cell has become a paradigm for 

^ understanding the molecular processes involved in gene expression 

and signaling within a cell. This system provides an example of a 
0^ genetic switch, as cells with identical DNA choose either of two cell 

cycles: a lysogenic cycle, in which the phage genome is incorporated 
into the host and copied by the host; or a lytic cycle, resulting in the 
death of the cell and a burst of viruses. The robustness of this switch is 
^ remarkable; although the first stages of the lysogenic and lytic cycles 

, ^ are identical, a lysogen virtually never spontaneously flips, and external 

^ stressors or instantaneous cell conditions are required to induce flipping. 

^ In particular, the cell fate decision can depend on the populations of two 

proteins, cl and Cro, as well as their oligomerization and subsequent 
binding affinities to three DNA sites. These processes in turn govern 
the rates at which RNAp transcribes the cl and Cro genes to produce 
more of their respective proteins. 

Although the biology in this case is well understood, the fundamental 
chemistry and physics underlying the bistability remains elusive. In this 
work, a dynamical model of the non-equilibrium statistical mechanics 
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is revisited, generalized, and explored. The low number of proteins 
and other sources of noise are non-negligble and corrections to the 
kinetics are essential to understanding the stability. To this end, general 
integral forms for advection-diffusion equations appropriate for finite 
element methods have been developed and numerically solved for 
a variety of mutants and assumptions about the state of the cells. 
These solutions quantify the probabilistic and flux landscapes of the 
ensembles' evolution in concentration space and are used to predict 
the populations of the cell states, entropy production, passage times, 
and potential barriers of wild type and mutant bacteria to illuminate 
some structure of the configuration space from which Nature naturally 
selects. 
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1 Introduction 



Life is an entanglement of lies to hide its basic mechanisms. 
William S. Burroughs, The Place of Dead Roads 

1.1 Computation, Noise and Limit Cycles in Biophysics 

It is a season in which the porousness of the borders delineating the subjects 
that comprise Knowledge are growing more so, or altogether vacuous. The 
overlap is revealing some truths to be quite universal. It has been said that more is 
different Anderson [1972] and information is physical Landauer [1961], and concepts 
involving emergence and entropy are proving to be fundamental to understanding 
the physical basis of biological systems. The notion of a complete and clockwork 
deterministic universe has been left behind, and the necessity of dealing with 
probabilities, partial information, measurement, and open systems has become 
essential to all of the major sciences. 

Among these disciplines, the common goals have become those of elucidating 
which states are possible, determining the relative probabilities or amplitudes of 
those states, and examining their stability and robustness. In the subject of physics, 
the need for these particular goals is most easily evidenced by the existence and 
necessity of quantum mechanics, which assigns to each system a wavefunction 
from a spectrum of permitted states; entanglement amongst these states; and, by 
measurement, results ultimately in the collapse of the wavefunction to a single 
eigenstate of the spectrum. In chemistry, it has become increasingly necessary to 
delineate the single pathway that a molecule may take, as it selects a seemingly 
random walk in a solvent, or geometrically aligns with a catalyst in a reaction. 
In biology, evolution provides a framework by which we understand how each 
member in a myriad of diverse species becomes finely adapted to their single, 
specific niche, and how these populations are interrelated. 

Although the underlying theories are generally accepted, there remain open 
questions and challenges in the implementation of these theories that are similar 
across disciplines. Analogies can even be present in the mathematical formalisms 
that are often developed in parallel for these subjects. For example consider time 
evolution in quantum mechanics, Brownian motion, and population dynamics. 
The integral of a specified initial condition is commonly implemented iteratively 
as a matrix operator acting on a vector which represents the quantum state in 
the case of quantum mechanics, the population densities in the case of chemical 
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dynamics, or the population density in the case of an ecosystem. 

The limit cycles are of considerable interest. Although teleological explanations 
in evolution are largely vanquished Mayr [2002], at least some of the properties of 
living matter can be predestined by the highly non-negligible physical and chemical 
constraints placed on it by the supporting planet. Understanding physically how 
systems acquire their properties through gene expression and feedback from an 
environment leads to fundamental physics. However large the number of unknown 
degrees of freedom, they serve as a source of energy and entropy for the system, 
or as a place to deposit energy and entropy. In some shape or form, this exchange 
introduces a mechanism for choice, chance, or randomness into the time evolution 
of the system, which is commonly implemented as a matrix acting on a vector.^ Be 
it the decoherence processes that select an eigenfunction out of entangled quantum 
states or the natural selection processes that select the fittest of a species, it is time 
that selects an actuality out of the sea of possibilities that could potentially occur. 
That is to say, whatever the equations, it is the asymptotics of those equations that 
are of interest. When the system is open, the nature of this selection speaks equally 
as much of the measurement device used or the planet underfoot as it can of the 
system measured or the species selected. 

Assuming that we have equations that describe these phenomena, is it feasible 
to imagine finding solutions to them? The question we ask is what the best models 
are for the physics we can develop given the computational resources that are 
available, and how we can improve the precision of the physics as our algorithms 
and computing power improve. When we point our computers at the abstractions 
we use for these systems, the polynomial algorithms seem nowhere to be found. 
While for biological problems, even by counting only the emergent degrees of 
freedom (proteins) of the emergent degrees of freedom (RNA) of the emergent 
degrees of freedom (DNA) of the emergent degrees of freedom (the nucleotides) 
of the emergent degrees of freedom (the periodic table) of the emergent degrees of 
freedom (the protons and the neutrons), our ability to directly probe the relevant 
scales of length and time is limited by our current information technology. 

Because many questions in one field are categorically identical to questions 
in another, the achievements and methods in one science are non-negligible to 
the others. Often they are the same achievement or require cooperation. The 
contribution of physics thus far to biology is largely that of engineering: it is the 
microscopes, spectrometers and transistors that are mainly of service. There is 
no evidence suggesting that the physics of biology requires more then quantum 
mechanics, statistical mechanics, and classical mechanics to operate, but the extent 
to which the solutions to our equations lag behind our equations provides ample 
room for doubt and countless theorems remain to be demonstrated to show this. 

^Consider, for example, x = Ax, ihdt^ = HY, dtP = —di{DijdjP — FfP) with solutions, x{t) = 

e^^x(O), Y(f) = e-f^ /o ™^Y(0), and P{t) = e^oi^^j^i-Pi)^^p{o) respectively. The essential problem, 
be it a continuous Markov process described by A, Schrodinger's equation, or the Fokker-Planck 
equation, the equation boils down to the evaluation of an exponential of an operator. 
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As an example, biology has provided us with the ingenuity through experiments 
to deduce the relevant underlying networks for the A phage Ptashne [1992]^, the 
flexibility of those networks Atsumi and Little [2006], Little et al. [1999], and how 
to introduce and parameterize the empirical force fields that are introduced to 
study the dynamics Shea and Ackers [1985]. These numbers are improved every 
day as better microarray data is accumulated and understood ^. 

Evolution, like gravity and thermodynamics, is so seamlessly entangled into 
our lives as to often go unnoticed. Once, however, the concept is understood it is 
quite self-evident, and is as readily seen as the selection of an item off of a menu 
or the selection of students or employees from a pool of applicants. This interplay 
between chance and determinism, the arbitrary and the necessary, and populations 
and their individuals provides more then enough room for complexity to emerge, 
even before granting a capacity for volition to the individuals. 

The key mechanism again is emergence - be it quarks to form a proton, nu- 
cleotides to form DNA, or cells to form an organism, the salient necessity for 
evolution to occur at any hierarchy of a universe is a mechanism to build up sys- 
tems of arbitrary complexity out of simpler and readily available resources from 
the nearby hierarchical levels Jacob [1977], McCoy [2010]. Unsurprisingly, we find 
ourselves present in such a universe, and Nature has no shortage of mechanisms 
to build complex systems, and clever algorithms for their improvement Darwin 
[2011], tol [2011]^. 

It has been said that the universe is a quantum computation Lloyd [2006], and 
it will not resist being perceived as one, but other modes of computation the 
universe readily performs, those using the counting of molecules, are worthy of 
deep reflection Bennett [1982]. These classical modes can also supply the speedup 
to which quantum computers esteem. This Brownian computer, as Bennett calls it, 
is not dissimilar to a quantum computer, in the sense that it relies on interference 
to decrease the probability of certain possibilities and increase the probabilities 
of others. The fundamental difference is in the nature of the interference, which 
for a quantum computation occurs at the amplitude level, where numbers can be 
complex, and for a Brownian computer this occurs at the probability level, where 
the numbers are strictly positive. The fundamental similarity between the two 
modes of computation is that they occur essentially in parallel. 

To illustrate, consider an impenetrable door, nay, there is an alphabet with four 
letters and one password in 4^ words that will let us in were we to enter that 
word. The keepers of the door also granted us a classical, quantum, and chemical 

^and references therein 

^In a just universe, the models such as discussed in this thesis would be coupled with data. The 
rates used in the partition functions should be fit simultaneously with the numbers used in the 
diffusion tensor as well. We do not take such a course here but provide the reference Kinney et al. 
[2007], in addition to Shea and Ackers [1985] to illustrate how AG numbers are accumulated 

^The notion of species-as-algorithm is evident in Toffoli [1982], where he gives the example of a 
typical control function: the octopus. 
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keypad. The classical keypad takes too long, qubits have their own agenda at those 
numbers, so how then will we ever open the door? 

What is a chemical keypad? Say we have n sites and we enter our code by 
binding 1 of our 4 molecules onto each site, and if we get the arrangement just 
right, then the door opens. A molecule will simply go away after some time if it 
did not fit. 

Our only hope then without any more prior information is BennettaAZs Brow- 
nian computer. We simply gas the keypad with concentrations of all the molecules 
and wait. Sooner or later the necessary configuration will occur and our key is 
found and the door is opened. 

This is still a classical computation in the sense that it is no different than our 
fingers punching numbers into a keypad as far as the information processing is 
concerned. Indeed the speedup here is simply of scale, as molecules can typically 
bind and unbind to sites and each other faster then we can push buttons. The 
probabilistic nature of this process comes only as a result of our description 
of the scenario. We have only macroscopic information, the thermodynamic 
variables, and may be able to tune them to some degree, but the goal is a particular 
conformation. This information is not accessible to thermodynamics proper, but 
we know the answer once the goal has been achieved, as we read it out of the 
keypad only after the door is opened. As is often the case, the key itself might fall 
out of the door and this information too is lost, but sometimes just an open door 
is sufficient to proceed. 

Evolution is such a mode of computation: a wide variety of possibilities are 
attempted, and if any open doors, or are built to last, then doors open and things 
last. The rate and mechanisms it achieves to do this provides curious insight into 
new computers and algorithms for achieving and dealing with computational 
complexity, typically involving non-binary modes of thinking. 

Granting a billion years from the birth of Earth to the birth of the first prokary- 
otes, say 4 * 10^ bp (base pairs) Mayr [2002], it can be infered that the effective 
rate at which natural selection processed the 4^*^^^ configurations is 4^*10^-27 y[z. 
This is a rephrasing of Levinthal's paradox in protein folding for the problem of 
evolution. It is thus safe to infer that natural selection is not picking numbers at 
random from a flat distribution. 

Genes are an example of an algorithm that has provided exponential speedup 
to this sampling. Genes might be a transcript for a protein and take up a thousand 
nucleotides of space, and at this stage Nature's trial and error processes of evolution 
revolve around a sampling of the space of gene combinations; this is a considerably 
truncated space from the space of permissible DNA chains. Sexual reproduction is 
one example exhibiting the benefits of such a sampling, and it emerged as a means 
of speeding up the quest for fitness about 2.7 billion years ago. The DNA of the 
resulting offspring is far removed in a symbol space from either parent's DNA, yet 
the genes persist. It is quite a luxury that a child should take after their parents. 

With a configuration space so vast and the demands of a niche so extreme, one 
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can hardly suspect a few random bits flipping to send one to a more favorable 
potential well. Quite the contrary is fact and indeed the structure of DNA, with 
its Watson-Crick pairing and strength of its bonds, is largely preventative of just 
this situation. Unsurprisingly, there exist proteins and mechanisms that specialize 
in DNA repair to further parry these crude attempts to debase the bases. In 
eukaryotes, DNA is further protected via a profound capacity for packing and 
elaborate protection from histones and proteins. 

So it is genes that drive development and the space of gene configurations that 
DNA prefers to play trial and error with, for exactly the reason to exponentially 
speedup the search for fitness. There is then a very strong and productive motive 
to treat the genes and their functions as the relevant degrees of freedom, and 
indeed that is the approach and success of systems biology. Ptashne and Gann 
[2002], Niisslein-Volhard [2006] 

It can be seen that evolution can proceed with large steps. There are mechanisms 
in which the step nature of evolution of DNA itself is evident. This revolves around 
the discussion of genes, not nucleotides, where very specific and relatively long 
transcripts can be inserted irreversibly into a genome. Such is the case we study 
now, the phage A infection of an E. coli. The biological content of this thesis is well 
known and actively extended, St. Pierre [2009], Ptashne [1992]. 

A is not alive and requires an E. coli host cell to provide the habitat and 
ingredients, in particular RNA polymerase and the ribosome, to express and 
replicate the genes present in that genome. A head provides protection to the 
DNA outside of a cell as the virus uses a tail to seek out a host upon which it 
pierces the cell membrane and injects the DNA into the cell. A special feature of 
the genome is the six bases on each end; they are unpaired 'sticky ends,' and are 
dual to each other. These sticky ends stick to to each other upon entering a cell 
and form a ring. The viral DNA then waits for the cell contents to regain their 
closest approach to equilibrium. This tricks the RNA polymerase and ribosome in 
the host into expression and synthesis 1.1 of its genes. 

This splitting of a cell population into different possible fates can be seen in a 
variety of circumstances. Some cells burst into a multiplicity of phages immediately 
(lysis), while the bacteria that do survive are irreversibly transformed into a new 
genome, now nearly 48,502 bp longer, and henceforth replicating along with its old 
identity and the entire phage (lysogeny), see 1.2. This behavior seems more akin to 
a jewel heist then a thermal, stochastic, or a quantum fluctuation. It is believed that 
noise influences the decision at the level of selection between fates. This noise has 
been studied elsewhere, Zhu et al. [2004], Wang et al. [2010], Arkin and Ross [1998], 
but the exact nature of this noise is not understood, and is typically presumed 
Gaussian for computational and analytic convenience. The methods developed 
to study this noise can easily implement more complicated noise, for example, as 
known functions of space and time. This switch is a perfect example of a switch 
component of a Brownian computer, but in the end just one component of a larger 
A gene network, presented in Arkin and Ross [1998], St. Pierre [2009]. 
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Figure 1.1: A phage Bacteriophage Borrows Machinery from the Host Cell 

requires a host for 



RNAp and ribo- 
some to express its 
genes. 




Turing ^s Dream Machines 



The rates at which these transcriptions occur enter the moders equation of 
motion as simply numbers and set the time scale of the computation, i.e. the 
expression of the gene, so let us revisit RNAp, and the ribosome, to first marvel at 
two of the most amazing molecules this side of the big bang, and chieftains among 
the hyperspecialized computational hardware inside a cell. 

Indeed it was RNA polymerase that was the source of inspiration for the 
Brownian computer in the first place, and its functioning had been described in 
information theoretic terms already in 1982 Bennett [1982]: 

To see how a molecular Brownian computer might work, we first 
consider a simpler apparatus: a Brownian tape-copying machine. Such 
an apparatus already exists in nature, in the form of RNA polymerase, 
the enzyme that synthesizes a complementary RNA copy of one or 
more genes of a DNA molecule. 

Bennett goes on to describe the RNA polymerase action in terms of a one- 
dimensional random walk along a DNA chain. In the forward reaction RNA 
polymerase opens the DNA to read the content and append to the growing RNA 
strand the nucleotide part of the relevant nucleotide pyrophosphate ATP, GTP, 
CTP, or UTP and subsequently releasing the pyrophosphate molecule into the cell, 
advancing a notch on the DNA chain. In the backward reaction RNA polymerase 
takes a free pyrophosphate (PP) from the cell and combines it with the last RNA 
bit and releases into the cell the ATP, GTP, CTP, or UTP. In equilibrium this reaction 
would balance and RNA and protein synthesis would not occur. However, the 
metabolic processes are such that a bias in the direction of RNA synthesis occurs 
by removing free PP and supplying ATP, GTP, UTP, and CTP to the vicinity of RNA 
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Figure 1.2: The two 
E. coli outcomes, 
lysogenic above 
and lytic below, of 
a A phage infec- 
tion. The Or is 
yellow, the active 
cl gene is blue, the 
active Cro gene 
is red, and grey 
are other active 
genes in that stage 
of the cell cycle. 
The Att gene in 
green is the gene 
responsible for 
incorporating the 
virus into the 
bacteria. White 
regions are DNA 
segments inactive 
during that stage. 

polymerase. The numbers used in his analysis are RNA polymerase advancing at 
30 nucleotides/s (330 forward reactions -300 backward reactions), and dissipates 
20kT per nucleotide all while making less than one mistake per ten thousand 
nucleotides. Bennett [1982] 

By 1999, this one-dimensional random walk was directly observed in single 
molecule experiments Guthold et al. [1999]. The rates therein were much lower 
(1.5 nucleotides/s ± 0.8 nucleotides/s), but these processes were likely slowed 
down by the experimental conditions which confined the DNA to diffuse along a 
2D substrate, restricting the motions of many degrees of freedom, to facilitate the 
observation of RNA polymerase moving along the DNA as if it were a ID chain. 
It had been argued that effectively diffusing in lower dimensions enables RNA 
polymerase to find the promoters much more easily. 

A has developed an additional strategy that speeds up this transcription in- 
volving cl protein Ptashne [2011]. It involves a relatively small cooperativity 
(-2kcal/mol) binding between a cl dimer and RNAp at the promoter site, which 
also shortens the time RNAp requires to diffuse to that point. This is a small 
number and can not alone describe the observed factor of 10-100 increase in tran- 
scription rate when cl is present. An explanation I provide here is that this speedup 
could be a result of another complex binding process, where cl octamers short the 
DNA loop, as has been observed. That this results in a speedup of transcription is 
a result of the simple fact that, if RNAp can maintain a 30 nucleotide/ s rate, this 
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rate is increased in units of m/s if the chain is shortened and there is simply less 
distance to diffuse to the end of the transcript and back to the promoter. Indeed the 
looping between the Or and the Oi leaves RNAp again approaching in physical 
space the promoter as it is finishing the transcription of cI. This increase does not 
occur for the Cro gene since the shorting of the loop does not contain the Cro gene, 
only the cl gene, and requires cl to be occupying the promoter site for Cro anyway 

The subsequent RNA strand is fed to the ribosome which implements a many- 
to-one dictionary, converting three consecutive nucleotides into an amino acid. 
The 64 4-letter words of length 3 get mapped to 20 amino acids, aiding in the 
formation of the protein. This allows for a large multiplicity of DNA words that 
can create the same protein and thus allow for a finer resolved selection process to 
occur amid the words that produce a given protein. 

cl is a chain of 236 amino acids and cro is 66. The rates at which these processes 
occur are actively measured and correspond to the value of v^i for cl and Cro in 
this situation. The rates for the lambda phage production of cl and Cro used in 
these simulations are 0.115 cl monomers /second for states with cl assisting RNAp, 
and 0.30 Cro monomers /second. Other impacts on the rates are the 0.01045 cl 
monomers /second contribuition from certain configurations, and that the ribosome 
produces 20 Cro for each reading of the Cro transcript. These numbers are chosen 
to correspond to the values use in Zhu et al. [2004], Wang et al. [2010] where this 
dynamical system is investigated. The complete set of parameters used in this 
simulation will be shown in the dynamics section. 

There are two sources of non-linearity in the equations of motion, first from 
there simply being multiple binding sites, but the main complication, and source 
of error in this work, is that both cl and Cro first dimerize in the cell and bind 
in dimer form to DNA with a helix-bend-helix motif. Different configurations 
have slightly different affinities and thus occur with different frequencies. All the 
configurations however have negative AG and as such are spontaneous reactions; 
they are energy favorable in comparison with an empty strand of DNA. In addition, 
cl dimers are found to bind cooperatively on adjacent sites. This added energy 
bonus from cooperation aids in the competition with Cro. Cro is also limited by 
the fact that both Ori and Or2 site need to be unoccupied for RNAp to bind to 
the Cro promotion site. The partition function is dominated by the triply liganded 
species and the faster rates and large bursts for Cro get largely repressed by the 
relative scarcity of the singly liganded states in comparison to the doubly liganded 
states that promote cl. 

The phage genome in a lysogen patiently waits for cell conditions to appear, 
a process known as induction, that switch the genes being expressed to phage- 
bearing ones and result in bacteria death and phage progeny. Ultraviolet light, 
small cell volume, and elevated Cro or depleted cl protein levels are known to be 
able to induce the switch. It is the object of this dissertation to further quantify the 
nature of this swith and compare with the words nearby in simple space. These 
are very specific mutations in the Or, chosen so that the binding properties can 
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be predicted, and switch properties studied. A result of this thesis, and a testable 
hypothesis, is that energetically the use of ultraviolet light for induction is a little 
extreme, as radiation targeted at the Or equivalent to 15kT should be sufficient 
to provide the energy for some of the lysogens to flip. Indeed, NMR techniques 
specialized to control cl to the extent that it can be removed from Ori and Or2 
should also be sufficient for induction even at larger concentrations. 

1.1.1 Chemical Oscillations 

That classical chemical reactions used to describe the balancing of morphogens 
within an embryo is sufficient to give rise to non-trivial spatial and temporal 
structures was noted carefully in the last works of Alan Turing Turing [1952]. See 
Kondo and Miura [2010] for a modern review. 

Classical treatment of chemical reactions is sufficient to understand many 
rhythms of the body, like a heart beating, and it is wise to push the utility of such 
treatments as far as possible. 

Reactions are not even necessary for oscillations to occur, indeed mixtures of 
perfect gases oscillate on grand scales in atmospheres in an attempt to maintain 
hydrostatic equilibrium for the mean density despite tidal or other forces Niemann 
et al. [2005], MAijUer-Wodarg et al. [2006]. 

The variety of chemical oscillations considered in this thesis are particular to 
those believed to be relevant in gene expression, which involve the binding of 
ligands to operator sites resulting in the creation of more ligands. 

The chemistry section of this work states quantitatively our basic assumptions: 
namely the partition functions used to describe the equilibrium statistical mechan- 
ics of the protein/ DNA reactions. It is stated to a reasonable generality what 
shapes these partition functions can make, facts that we will need to motivate 
the extent of accessible dynamics to which these binding polynomials along with 
transcription factors can give rise. 

At low numbers, chemical oscillations take on much greater complexity, where 
noise can play dominating roles. Most importantly this work relies on the discus- 
sions in Wang et al. [2008]. 

2.2.2 Dynamics 

Universal computation takes on a whole new meaning in chemical and quantum 
computation. One need not seek only classical answers to classical questions, 
as analog circuitry and quantum computation can greatly enrich the classes of 
problems that can be understood. The dynamics resulting from gene expression 
could be of interest to the purest of mathematicians, and in some sense are the 
n-body generalization of Hilbert's 16th problem, which remains unsolved. We 
applaud chemistry and biology for their courage in asking such questions. Hence 
it is quite appropriate to begin with the most general equations for dynamical 
systems. The qualitative and quantitative theories of dynamical systems are old 
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and advanced Liapounoff [1893], Volterra [1931], Arnold [1978], Perko [2001]^ and 
the need for general methods to general equations in simulating gene expressions 
is noted in Zhu et al. [2004], Wang et al. [2010], Aurell and Sneppen [2002]. 

The codes used for this work were written to maintain this generality and the 
software can be used to study dynamical systems over real or complex fields. 

1.1.3 Diffusion 

Of the few equations Schrodinger reluctantly allows appear in What is Life? 
Schrodinger [1992] the first is 

dfP = DAP (1.1) 

Which he need not trouble you to explain. The gradients and diffusion of which he 
speaks are in physical space, whereas we will be concentrating on concentration 
space. The motive and solutions are, however, similar. 

The robustness of the phage A is of the type Schrodinger would have claimed 
at the time to be impossible. He deliberates extensively on the ^/~{n) law and 
expected that where there is order, such as ferromagnetic order, there will be a 
requirement of a large number of particles involved. Yet the order one would 
expect from large concentrations is still present at the small numbers of 100-300 
proteins and has given rise to the so-called stability puzzle Aurell et al. [2002]. 

The analogous equation of (1.1) for the case of concentration space, to which 
numerical solutions herein are devoted to address, is 

dtP = di{DijdjP-FiP). (1.2) 

The codes allow for F/ to take on the generality prescribed in the dynamics section 
and Dij can be functions of space and time. Our mathematical analysis will proceed 
with the simplest case of a positive and diagonal diffusion tensor. 

A main contribution of this thesis is the adaptation of (1.2) for use with finite 
elements. 

2.2 Computers and Phage A 

Before we posit equations and automate solutions let us note that biological 
systems are not dissimilar from computers and their behavior is not dissimilar 
from a program. There is no hope of putting it more eloquently then Rothstein 
[1982]: 

aAYComputers and living systems share many characteristics: they 
behave selectively, show tremendous differences in response to similar 

study of dynamical systems, to which noise is largely a perturbation, can on no accounts be 
trivialized or neglected. Many texts on the fundamentals are available, in particular Perko [2001] 
illustrates the linearization methods to understand the local, qualitative dynamics of a system. 
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inputs, and are nonequilibrium systems, generally metastable. They 
usually require an ongoing dissipation to maintain their characteristic 
behaviors, which depend on stored information. Both involve loosely 
coupled subsystems, use elaborate internal communication and control 
systems, and key interactions within them and with their environments 
are typically irreversible and all or nothing (nonlinear)/ 

One may also add that they have already been selected from a vast space. 

The phage A is simply one of the 4^^^^^ molecules DNA could have chosen 
from this vast space. Yet, this sequence warrants its own name, for it is uniquely 
adaptated to secure the fidelity of that same sequence through the infection of a 
host cell. To do so, it chooses selectively between lysis and lysogeny and requires 
an ongoing interaction with RNA polymerase and the ribosome to achieve its 
steady state protein concentrations through the transcription of the cl and Cro 
genes. This involves the coupled subsystems, the cl and Cro proteins and the 
Or, with elaborate internal communication through oligomerization processes, 
and key interactions within them are indeed nonlinear. The fate of the decision is 
thermodynamically irreversible. Let us elaborate. 

1.3 phage A 

The enterobacteria phage A infection of Escherichia Coli is proving to be an insight- 
ful example in biology of how the small-scale molecular interactions involved in 
gene signaling and expression lead to distinct outcomes Ptashne [2011], Ptashne 
and Gann [2002]. In this case, these outcomes are the very life and death of the 
cell. A history of ingenious experiments has revealed the interactions involved in 
this determination of cell fate in exacting detail Ptashne [1992], and the A phage's 
strategy for reproduction is a clear example where the molecular strategy involved 
in the evolution and natural selection of the organism are well understood. The 
questions then remaining are the chemical and physical questions in the cellular 
environment, particularly, what the interactions between these macromolecules 
are, and the type of model that is best used to quantify these interactions. The 
timescales involved in the lambda phage genetic switch are on the order of 10^ 
s to 10^ s, thereby placing the system beyond the limits of an all-atom molecular 
simulation with an empirical forcefield. Phenomenology has been introduced, and 
the origin of the model herein is the formulation of Shea and Ackers [1985] with 
subsequent revision in Aurell et al. [2002] and propagated by Wang et al. [2010], 
Ao [2004] and in the framework of Wang et al. [2008]: the diffusion equation. This 
has been adapted herein for the finite element method. 

The first description of the A phage system as a two-state switch was given by 
the biologists Lwoff, Jacob, and Monod Lwoff [1972]. The A phage is a virus that 
recognizes and binds to its host, the bacterium Escherichia Coli, and subsequently 
injects its genetic material into the bacterium in the form of double-stranded DNA, 



Chapter 1. Introduction 



13 



forming a ring. After the virus has infected the cell, there are two main mechanisms 
by which the viral DNA may be replicated: the lytic cycle, or the lysogenic cycle. 
In the lytic cycle, the bacterial DNA is destroyed, and the virus takes control of the 
metabolic activity of the cell with the goal of replicating its own DNA. The host 
cell is quickly overcrowded with new viruses, and the cell bursts, releasing the 
new progeny viruses. In the lysogenic cycle, the viral DNA is incorporated into 
the DNA of the host, and cellular replication occurs normally. The host cell is not 
destroyed, and the virus does not produce progeny. A summary of the genome is 
depicted in Fig. 1.3. The goal of the experiment of Lwoff et al. was to understand 
the genesis and predictability of the lytic cycle. In their work, a collection of E. coli 
cells infected with the A phage was irradiated by a moderate dose of ultraviolet 
light. In contrast to the expected scenario observed without radiation, in which 
many generations of cells underwent uninterrupted reproduction, the irradiated 
cells stopped growing and 90 minutes later a large fraction of the cells burst into a 
crop of A viruses. Repetition of the experiment led to reproducibility of the result, 
and the phenomenon was understood to be an example of the switching on or off 
of a gene controlling the lytic or lysogenic fate of the cell. 

Figure 1.3: Some 

C'l o Ti'iiiisciipt 

basic features 



of the A phage 
genome. 




The next major step was the discovery of the Or region Ptashne [1992]. The 
48502 base pairs of the genome had been sequenced, and experiments showed that 
between basepairs 37940 to 38020 there existed a region, the Or, between the cl 
and Cro genes where the cl and Cro proteins compete on binding sites for their 
transcription. This region overlaps with two RNA polymerase promoters, Prm/ Pr/ 
for promoting cl and Cro, respectively. Prm overlaps with Or^ and Or2, while Pr 
overlaps with Ori, Fig. 1.4. This has given the cl protein the name A repressor, 
since in lysogeny the cl protein actively represses the expression of the Cro gene 
by blocking the Pr binding site from RNA polymerase. When Cro is present the 
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affinity is hightest to the Prm site, blocking cl promotion. 




The Or Region is a Genetic Syvitcli 



Figure 1.4: Zoom- 
ing into the Or site 
reveals more struc- 
ture, the Pr site 
overlaps with Ori 
and Or2 and Prm 
overlaps with Ors. 
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In this study, we focus on the fate of the E. coli bacterial cell at the moment that 
infection with the A phage occurs. 

In the chemistry section of this work, we introduce the equilibrium assumptions 
used in the model and tabulate the rates involved in some of the system kinetics, 
namely the binding of dimers to DNA, and the dimerization of the proteins in the 
cell. The AG values are the result of fits to experimental studies and are largely 
the focus of experimental work on which this model depends. All numbers have 
been chosen to correspond with the numbers implicit to Wang et al. [2010] which 
we repeat here. With modern crystallographic studies of these interactions and the 
state of empirical force fields, these values could in principle begin to be given a 
basis without the biological experiments. This avenue will not be pursued here, 
alhough it is representative of the type of hurdle involved in getting clinically 
useful information out of molecular structure-based studies and we stress the 
importance of such calculations to provide some atomic and molecular justification 
of this phenomenology. 

In the dynamics section of this work, we review the phenomenology that has 
been introduced. This phenomenology is interesting in its own right as a method 
for exploring non-equilibrium statistical mechanics behavior over large time scales 
given some understanding of the equilibrium picture. The success with the A 
phage is suggestive of rather general and powerful methods for exploring other 
gene regulatory circuits and complex systems as well. The essential content of the 
dynamical system of Shea and Ackers [1985] is contained in the rate equation that 
follows 



(1.3) 
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Figure 1.5: De- 
picted are three 
possible states of 
the Ori site. Each 
of the three DNA 
sites can be either 
empty, occupied 
by a cl dimer, or 
occupied by a Cro 
dimer. 



States of th e Site 





Here, the time derivative of the Xf protein is written in terms of the rate, v^., at 
which RNAp and ribosome can actively create more X/. Zx- is itself the sum over 
states that allow a creation process and is weighted by a sum over all states. In 
the cell the protein decays at a rate of . . No summation convention is employed 
with this equation. 

The inevitable limitation of Eq. (1.3) is the presence of noise, and corrections 
are necessary. To this end, the diffusion section of this work introduces an integral 
method for solving a class of partial differential equations including the Fokker- 
Planck equation dtP = —di{FiP — DjjdjP). This is to provide a numerical method 
of solution to the general formalism for using diffusion equations to describe noise 
in chemical equations, Wang et al. [2008]. This is an analogy to stochastic methods 
for chemical equations Ge and Qian [2010], Qian et al. [2002], and path integral 
methods for chemical equations Wang et al. [2010]. The class is then restricted to 
the case of flows induced by the dynamic assumptions of the phage. This method 
is typical of a finite element method and its origins even predate the existence of 
the computer that has made it so practical Logg et al. [2012]. The following two 
equations are the main result of this thesis. 



(1.4) 
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L = J^i^ut - —div[DijdjUt - FiUt])dx. (1.5) 

The solutions Ut and / = DjjdjUt — FjUt are the probabilistic and flux landscapes, 
respectively, that appear in the title of this work. 

a = L is a linear matrix equation in a function space when one specifies an initial 
condition. The dimension of the space depends on the computational resources 
at our disposal and the flexibility in this direction is just one luxury of FEniCS 
software.^ The expression a is bilinear, in Uf^st ^rid in v, when used for propagating 
time in the positive direction. The right hand side in this equation is only linear, 
depends only on v, i.e. Ut is fixed by the previous time step and thus by the initial 
condition. 

This is done by casting all the functions into a function space, or vector function 
space as needed. When these expressions hold for all v in the function space, Uf^^f 
is a solution to the related Fokker-Planck equation if Ut is. Integral nth passage time 
equations are similarly introduced and first passage time equations are solved for 
all the mutants. The integral equation corresponding to F/3/fi(x) — D3/3/fi(x) = 
— 1 follows 

/ {t{x){Fiditi{x))-Ddit{x)diti{x))dx = - [ t{x)dx, (1.6) 

/ {t{x){Fiditn{x)) — Ddit{x)ditn{x))dx = —n / t{x)Tn-idx. (1.7) 

The passage times proceed recursively. After the first passage times are found 
and inserted for Tq the second equation is rendered linear and the higher moments 
follow similarly. The situation in the Or requires only the first passage time; the 
extreme value of these numbers compared to the cell cycle is how the robustness 
is quantified here. There is some dispersion in the literature as to what that value 
is, as it is difficult to measure directly, but an estimate of the upper bound is that 1 
in 10^^^ flip in a cell cycle, see Little et al. [1999]. 

The integral equations are discussed in the diffusion section and further tab- 
ulation of the results occurs in chapter 6. The methods were kept general and 
modular to quickly facilitate deeper exploration as more computational resources 
are employed or parameters changed. It is taken for granted that the models of 
systems that evolve should likewise undergo evolution. 

The A phage also meets the criterion for a selective system to be worthy 
of intensive study; the A phage is a modular system, a system of interacting 
components that can be systematically permuted Rothstein [1982]. 

He also motivates purely computational approaches to theoretical biology. 

aAYProgress in theoretical biology, we believe, needs help from 



^Various mesh sizes were used in this thesis, a 600 by 600 space of quadratic elements, corre- 
sponding to a vector space of ^ 720000 dimensions was about the limit of a desktop computer. 
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more than physics and biology. It needs help from computer science, 
where any [selective system] that can be described by either a statistical, 
mathematical, or logical theory can be modeled. In particular, reconfig- 
urable systems consisting of interacting subsystems should be studied 
intensively, a AZ 

The A phage has proven to be a case in point of the benefits of such intensive 
study. The effects of some mutations in the Or have already been experimentally 
realized. Little et al. [1999], and A121, A323 and A323/ have been experimentally 
synthesised and studied. The most successful exploration in this varietal has 
been by far directly from the biologists, where DNA and proteins are directly 
manipulated and observed. It is the authors hope that in time the contributions 
from physics, mathematics, computer science, and chemistry will provide the 
means to observe much more quickly and cost effectively some of the effects 
possible in this vast configuration space, and direct biological experiments towards 
interesting possibilities. It is largely in the direction of this ideal that this thesis 
explores, and we will simulate those and similar mutations, proposing theoretical 
predictions to feasible experiments which could further constrain the models. 

lA Evolution 

Evolution, like gravity and thermodynamics, is so seamlessly entangled into our 
lives as to often go unnoticed. Once, however, the concept is understood it is quite 
self-evident, and is as readily seen as the selection of an item off of a menu or 
the selection of students or employees from a pool of applicants. This interplay 
between chance and determinism, the arbitrary and the necessary, and populations 
and their individuals provides more then enough room for complexity to emerge, 
even before granting a capacity for volition to the individuals. Much confusion 
about evolution can be avoided by clearly seperating what one is selecting from 
and selecting for, Mayr [2002]. This thesis takes leave from the noblesse oblige, in 
speaking of matters in which one is not an expert,^ but the motivation of such an 
investigation as this is largely rooted in Evolution, and such a discovery warrants 
discussion in a physics thesis. Indeed in Ao [2008] it is speculated that Evolution is 
the theory of everything physicists have vainly sought after. In addition, the main 
methodology of this thesis is one of tinkering; to do lots of trials and make a lot 
of errors. This process is in correspondence with how natural selection naturally 
selects, see Jacob [1977]. 

1.4.1 Natural Selection from... 

We will be considering natural selection from the set of DNA configurations. We 
choose 27 of these configurations corresponding to mutations in the Or region 

^somebody has to bridge the gaps...Schrodinger [1992] 
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for which present knowledge enables a prediction of what the binding affinities 
would be. It is demonstrated that just this restricted space is enough to sample 
both sides of a bifurcation. 

2.4.2 Natural Selection for... 

If biological systems are like other physical systems they seek to minimize energy 
and maximize entropy. We will rank our mutants from fittest to least fit based on 
passage times, energy, and entropy. Natural selection is a fine balancing of many 
considerations. 

1.5 Three modes of Computation 

Many modes of computing with molecules have been proposed and these can be 
split into three overlapping categories. 

Although the number of genes in the human genome has likely been constant 
for a few million years, the number in the refereed literature steadily decreases 
with time. In the human genome there are believed to be around 25000 protein 
coding genes, and much more of the genome is understood to be dedicated to 
regulation and signaling processes than initially thought. 

Cooperativity between regulatory proteins and RNAp for transcription is 
seeming to be the rule rather then the exception Ptashne [2011]. The signaling 
involved masterfully involves the time domain and is the major port of departure 
between DNA type computation and a Turing tape. Whereas a Turing maching 
requires a head to physically transform the bit strings along the tape, DNA 
strictly avoids this phenomenon and instead prefers to modify the tape through 
binding of ligands. This is the tragic flaw of Bennet's proposed DNA computer in 
1982 in which a ribosome like structure alters base pairs. It may be more easily 
implemented on RNA then DNA, but still the lesson to be learned is to program 
with the time domain as much as possible. The paradigms of causality are being 
replaced by the paradigms of synchronicity; GPU, DuoCore, QuadCore, quantum 
computation, and chemical computation are all parallel computations that rely on 
a balancing of rates in order for theircomponents to sync and perform properly. 

1.5.1 Binding Sites 

Anything that can count can compute, and RNAp transcribing a gene into a 
protein is such a means of computing. This is the type of computation discussed 
quantitatively in this thesis, using the phage A as a concrete example. In the 
dynamics section we will discuss the rich space of analog computation, which 
computation with binding sites can simulate. 
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2.5.2 Polymerization 

A type of computation DNA can employ, in which a large configuration space is 
rapidly sampled through polymerization processes has demonstrated working, 
fast algorithm for NP problems. Adleman [1994, 1998], Braich et al. [2002] 



1.5.3 Solitons 

A method of computation that macromolecules could employ is through solitons. 
The effective field theories for some of the empirical force fields used to describe 
the molecular dynamics are known to emit soliton solutions, Yakushevich et al. 
[2002]. The presence of soliton solutions is sufficient to define a bit, and the 
propagation of the soliton is a sufficient means to transfer a signal. The empirical 
model used in protein folding is similar in structure to the empirical models used 
in superconducting circuits. Cascading logic implemented through solitons in 
superconducting circuits is the most promising approach to superconducting- 
based computation, Likharev and Semenov [1991], Semenov et al. [2003], Ren et al. 
[2009], Ren [2011]. It will be found that biological systems have the same physical 
limitations as engines or other computational devices and deep understanding of 
those limitations are required, Likharev [1982]. 

To specify, consider the Hamiltonian used in Yakushevich et al. [2002], 

TT V^/^T J.2 , 1 T i2 , • 2/^n+l,l ~ 

^ = ^M^n,l + o ^n,2^/>n,2 + ^n^Sin^{^ ^ ) 

^ (1.8) 

+ en,2Sin^{- \ ) + yocB{<Pn,\><Pn,2))- 



Now consider a Hamiltonian used in superconducting circuits Averin et al. 
[2006], 

^ = + 2E;Sm2(|i) + El(<^„+i -cpn- cp^^f]. (1.9) 

We can see that the correspondence is not just a matter of romance. Letting the 
nucleotides correspond to junctions, and the phase of the superconducting order 
parameter correspond to physical angles, one can proceed quite formally. 

The solitons would then be of the sine-Gordon varietal and approachable 
through standard methods, Rajaraman [1987]. 

It is a mystery of gene expression and cell behaviour how genes turn on and 
off, how several distinct and large regions of DNA, genes, can be quiet for the 
majority of a cell cycle and then be triggered for expression of certain proteins 
at critical moments. Although these processes are often known to be caused by 
chemical triggers, mechanical triggers are not ruled out, and anything that can 
induce long-range spatial and temporal correlations over a chain of DNA should 
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be considered deeply. Solitons typify long-range spatial coordination, and are at 
least an entertaining possibility of signaling modes available in DNA. 

The analogy is suggestive also of modeling systems via direct mimicry Junc- 
tions could be chosen with capacitance and junction energies to directly simulate 
on a transmission line the dynamics of a DNA chain. Similar circuits can be 
designed for protein folding with a corresponce between the junctions and the 
amino acids. Although DNA and proteins are brilliant, they are not particularly 
fast with that brilliance, and such circuits could boast protein folding faster than 
proteins can fold a protein. The lesson to be learned in superconducting circuits 
from DNA is a promise of diverse effects emerging from the aperiodicity of a 
crystal. One might be capable of developing dramatically sophisticated circuits, 
without a large number of components, by using similar-but-different junctions. 
The richness of the protein functions themselves are believed to be derivative of 
their structure, which itself is yet another chain, this time of geometric motifs, 
helices, hairpins. One is left curious to find analogues of these structures in the 
circuits. 

1 . 6 Methods of Numerical Analysis 

The rich diversity intrinsic to biological systems lends necessity towards mathe- 
matics and computation that is equally rich and diverse. In computational terms, 
biophysics necessitates modularity, as code needs to be useful and applicable 
in a wide variety of settings. Emphasis has been placed on exactly this point 
throughout the design of numerical investigations reported in this thesis. This 
outlook has allowed for scanning large regimes of parameters of particular systems 
as well as investigating many systems and models. It is a well known problem of 
computer science that generality comes at the cost of requiring large amounts of 
code to consider all possible cases. This problem has largely been circumvented 
through the use of automated form compilation, in which the code written is 
designed to write more specific code when a particular instance is requested, Logg 
et al. [2012]^. Such a method lends to an elegant correspondence between the 
mathematics describing a problem and the corresponding computer instruction. 

The mathematics employed for this analysis is a modern incarnation of finite 
element methods developed before computers existed to provide a general method 
to solve or approximate large classes of partial differential equations from theories 
of hydrodynamics, electromagnetism, and elasticity, and is sufficient to solve the 
partial differential equations arising from biophysics as well. We find that even 
in circumstances where analytic solutions exist and can be found, numerics are 
more illuminating then the resulting expression and provide greater accesibility 
to the particular values of a function at given points in the domain. An analytic 



^One cannot overstate the importance of the software and contributors therein. Their achieve- 
ment is momentous. 
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expression may require just as much numerical work to, say, evaluate and sum the 
necessary hypergeometric functions when a particular value at a particular point 
for a particular set of parameters. Thus, even in instances where symbolic analysis 
of solutions are accesible we insist on numerics. 

Another important necessity of computational techniques used to study bio- 
logical systems is scalability. Computer instructions, if not the computer itself, 
should easily adapt as numbers of parameters, degrees of freedom, and resolution 
of meshes increase. Equally as diverse as the biological systems themselves are the 
assertions found in the literature about the most relevant degrees of the system, the 
rate of changes exhibited by the system, the parameter values the system exhibits, 
etc. Scalability ensures that the choices made in approximating a system can be 
based on physics or other considerations and not depend exclusively on access 
to computational power. For example, better computers should result in better 
approximations but not be a prerequisite for the utility of the method; indeed, 
the use of these methods predates the existence of the computer. The desire for 
scalability also motivated the use of FEniCS software which can be parallelized for 
larger systems and meshes. The scaling that these methods utilize are suboptimal 
in that as the number of system degrees of freedom are increased, the computa- 
tional resources required for the same resolution increases exponentially. For an 
optiplex duo core, a computer of modest price and performance, we found this to 
limit the size of mesh to two dimensions with 600 quadratic Lagrange elements 
in each direction. Three-dimensional queries required the use of adaptive mesh 
techniques to locate the active regions of interest and localize the resolution to 
these areas. The computer instructions, the human interface to the computer, scale 
ideally. That is to say, the only instructions to change as the systems analyzed 
change are the definitions of the system. 

This partial scaling manifests in the code simply crashing when more resources 
are requested then available, but would continue to work on machines or clusters 
where the resources are indeed available. A pleasant corrolary of the active 
development of FEniCS is that code already written ages like a fine wine. The 
benchmarks stated have changed over time, constantly speeding up, and one 
only needs patience for the physical interface to be as fast as one desires, as the 
human interface is nearly optimized. If you want a solution, it is now as easy as a 
well-defined question. 

Another important aspect of the computation, to the author, is that all results 
and methods are reproducible and extendable from open-source utilities. This 
facilitates the transfer of knowledge from the more fortunate to the less fortunate 
and simply evolves faster than closed source analogues. That is not to say that 
open-source is better than closed-source so much as that it will be better, and con- 
tributions in the direction of automated form compilation are quickly eliminating 
the gap between the quality of freely available symbolic mathematics code and the 
privatized software for the same tasks. 
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1.6.1 FEniCS 

The diffusion equations herein would never have been solved without the expertise 
and their freely available software. If the codes used throughout this thesis are to 
be used succesfuUy a basic understanding of not only object-oriented programming 
is necessary, but also understanding and installation of the software available at 
www.fenicsproject.org is necessary. 

Data aequatione, spoketh Newton, quotcunque fluentes quantitae involvente flux- 
iones invenire et vice versa. Which Arnold put into modern mathematical prose: It 
is useful to solve differential equations. Finite elements do just that. 

1.6.2 ROOT 

ROOT is an open source utility developed at CERN and the standard analysis 
framework used at all major particle accelerators of the world, CER [2011]. The 
object-oriented tools like graphing and histogramming utilities largely extend the 
utility and accesibility of data generated in this thesis. The tools long since built 
for analysing nuclear data are readily available for the analysis of other data as 
well and provide automatable procedures for problems such as fitting peaks and 
simple accounting. 

In particular, a pyROOT paradigm was embraced to combine this functionality 
with the ease of python and integrate seamlessly with the python tools built by 
FEniCS. 

2.6.3 Aleph 

It is taken for granted that models of systems that evolve should evolve as well, 
and Launchpad documents the code and the evolution of this code, Borggren 
[2011]. 

1.7 Overview 

This work is largely a survey; it is the study of a biological system with the 
techniques of applied mathematics, engineering, computer science, and chemistry, 
directed at an audience of physicists and philosophers. As such it runs the risk 
of being neither here nor there, to shoot three bears while only infuriating them. 
To this end offense is inevitable. An ISI Web of Knowledge query reveals 170,303 
results for Finite Elements, 11,445 for A Phage, and 3,151,834 for genes. To contrast. 
Quantum Mechanics appears 31,673 times and Quantum Hall Effect 6,557 times. 
Thus, the references are largely a matter of taste, and do not attempt a reflection of 
history, so much as what I felt to be instructive accounts of the systems and useful 
ports of entry to the historical literature. 
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A goal is maintained to provide a perspective on generality while keeping the 
A phage as a concrete example. This work seeks to explain to a physical audience 
the details of a biological system. 

The philosophy of physics and information took inspiration from many places: 
Anderson [1972, 2000], Bennett [1982], Toffoh [1982], Likharev [1996, 1982], Lan- 
dauer [1961], and Rothstein [1982]. 

In matters of solitons there is Yakushevich et al. [2002], Averin et al. [2006], 
Likharev and Semenov [1991], Likharev [1996]. 

For background and data in the A phage, see Ptashne [2011, 1992], Ptashne and 
Gann [2002], Little et al. [1999], Atsumi and Little [2006]. 

For evolution there are good introductions, see Mayr [2002], Jacob [1977]. 

For mathematics, Arnold [1978], Perko [2001]. 

For early modeling and the head node for the taxonomy tree of the AG values. 
Ackers et al. [1982], Shea and Ackers [1985]. 

Later modeling, Aurell and Sneppen [2002], Aurell et al. [2002], Zhu et al. [2004], 
Wang et al. [2010]. 

The methods and ideas herein are largely in correspondence to: Wang et al. 
[2008], Zhu et al. [2004], Wang et al. [2008, 2010], Ge and Qian [2010], Feng and 
Wang [2011] and have benefited from numerous conversations with Jin Wang. 

An accomplishment of this thesis is one of method; it is the computational 
framework in Borggren [2011], this is but a modest interface and synthesis of tools 
Logg et al. [2012], CER [2011]. 
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Binding Polynomials 
of the Oji Region 



and the Chemistry 



Pray he under no constraint in this house. This is Liberty-hall, gentlemen. 
You may do just as you please here. 

Oliver Goldsmith, She Stoops to Conquer 

11 AG Measurements 

In this section we consider the partition functions used to describe the chemical 
binding interactions in the region that must be posited before we proceed to 
introduce dynamics. There is some choice so we explore multiple possibilities. 

For one, the original partition functions include explicitly RNA polymerase at 
a fixed concentration, whereas many new pictures also include bindings to the 
Ol region. A complication is that both cl and Cro proteins bind to any of the 
3 sites, on or Ol in their dimer form. It is also now established that up to 
octamers can be present and, amazingly, cl octamers can bind at once to the 
as the Ol and short the DNA loop exactly across the cl gene. We will continue to 
use the equilibrium relationship between dimer and monomer concentrations in 
our partition function, Reinitz and Vaisnys [1990], Wang et al. [2010], Aurell and 
Sneppen [2002], Zhu et al. [2004]. 

Again, x, y are cl, Cro protein numbers and [cI],[Cro] denote the dimer con- 
centrations. The simplest partition function used does not include RNAp binding 
affinities or concentration, and consists of the 27 states permitted by having noth- 
ing, cl, or Cro on each of the three sites. Of these states, 1 is zero liganded, 6 are 
singly liganded, 12 are doubly liganded, and 8 are triply liganded ^ The nature of 
the partition function is perhaps most evident when written in the form 

^ Since AG values are negative and the individual processes occur at similar rates to each other, 
this establishes a sort of band structure; the triply liganded configurations occur the most often, 
then the doubly liganded, then the singly liganded, and then the scarcest event is to find an open 
DNA configuration. 
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: = 1+ 

[cl]{e ^ +e ^ +e ^) + 

1/ ^^002 ^^020 ^^200 . 

[Cro]{e RT +e Rt + e ) + 



[clY{e RT +e RT +e Rt ) + 

i9/ ^^022 ^^202 ^^220 . 

[CroY{e RT +e rt + e Rt ) + 

[cl] [Cro] {e R^ + e ^ + e ^ + (2.1) 

^^201 ^^210 _^^120x 

e RT +e RT +e Rt ) + 
[clf{e-^) + 
[Crof{e-^) + 

[cJ]^[Cro](e RT +e RT +e Rt ) + 

[cJ][Cro]^(e RT +e RT +e rt ) 

This form makes it clear that it is symmetric with respect to permutation of the 
sites. The argument of Aurell and Sneppen [2002] is that RNAp needs an open 
site to promote Cro and an open Prm site to promote cl, and if the site is open, 
RNAp will be overwhelmingly first in line to bind and promote such states. It 
is known that the Pr site overlaps with Or-^ and Or^ and Prm overlaps with Or^. 
This symmetry is not present for the subspaces; the sum over the substates with a 
free Prm can promote x, a cl monomer, and the sum over the substates with a free 
Pr site can promote y for a Cro monomer 



1+ 



AG, 



010 



[cl] (e RT +e 



AG 



100 . 



RT 



) + 



AG, 



1 / - - 020 

[Cro] {e RT +e 



AG' 



200 . 
RT ) + 



[Cro]^ 



AG 



110 



RT -\- 

_ AG220 
RT -\- 



1 + 



r Tir^ 1/ ^^210 _^5l20x 

[c7][Cro](e RT +e rt ), 



AG, 



r -ri -—001 

[cl\e RT + 

T ^^002 

[Cro]e 



(2.2) 



with 
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[cl] = [x*pnc]/2 + e^^^^/^^V8 

^{AGci/RT) ^2AGci/RT (2.3) 

- ([^ * P^c] g + ^4 )^ 

[Cro] = [y*pnc]/2 + ^^^cro/i^T/g 

^(AGcro/RT) ^2AGcro/RT (2.4) 

-([y*p^^] — g — + — ^4 — ) 

Interlude: On dimerization.^ 

The key feature to notice here is that, since all the binding affinities are nega- 
tive, the sites are dominated by the triply liganded terms. The occurrence of an x 
promoting state is already rare, and a y promotion state exponentially more scarce 
then those. The partition function here is seen to be a sum over monomials of 
the forms, 1, [cJ], [Cro], [cJ]^, [Cro]^, [cI][CroUcI]^,[Cro]^,[cI]^[Crol and [Cro][cI]^. 
In the original Shea and Ackers contribution we note that additional monomi- 
als of the form [RNAp], [RNAp][cI], [RNAp][Crol [RNAP][cI]^, [RN AP][Cro]^ , 
[RN AP][cI][Cro] and [RNAp]^ were also included. This contributed to 13 more 
states included in the partition function for a total of 40. Only 40 of the 3^ states are 
permissible since an RNAp at Pr precludes any binding to Or^ and Or^. The — AGs 

2 

On Dimerization 

1 + 1 = 1 

Una goccia piu una goccia, fanno una goccia piu grande, non duel 
the Tarkovsky equation and explanation, Nostalghia 

Todos losfuegos elfuego 
Julio Cortazar, Todos los fuegos el fuego 

The dimerization equation used in this work originates from a quasi-steady state assumption 
for the dimerization reaction. 



[ckf = [cl2]e^' Wci] = [ch] + 2[ch] (2.5) 

This gives a quadratic equation with roots giving the cl dimer concentration, [CJ2] in terms 
of the total protein concentration, [N^/]. 



The root with the positive sign is the expression stated in Wang et al. [2010], Ao [2004] and 
the negative sign corresponds to the expression used in the calculations of Reinitz and Vaisnys 
[1990], Ao [2004], Wang et al. [2010]. The negative sign seems to be the better physical choice 
in that, in the limit of zero monomer concentration, one also has zero dimer concentration. 
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Ai23, wild — type Or affinities (kcal/mol) 


AG,-yfc 


protein 


Ori 


Or2 


Or3 


ref. 


equation 


monomial 


premotes 


AGooi/AGoio/AGioo 


-15 


-13 


-12 


Zhu et al. [2004] 


AGooo = 0.0 


1 




AG002/AG020/AG200 


-18.4 


-17.1 


-19.5 


Zhu et al. [2004] 


AGooi = -12.0 


[cl] 






AGoio = -13.0 


[cl] 




AGcoop = -6.9 


AGioo = -15.0 


[cl] 


Zx 


pnc = 1.5* 10~^^ 


AG002 = -19.5 


[Cro] 


Zy 


mpi = 3^ 


AG020 = -17.1 


[Cro] 


Zx 




AG200 = -18.4 


[Cro] 


Zx 


AGij]^{cont) 


AGoii = AGqio + AGqoi + AGcoop 


[cl]' 






AGiio = AGqio + AGioo + ^Gcoop 


[cl]' 


Zx 




AGioi = AGioo + AGqoi 


[cl]' 

L J 






AG022 = AG020 + AG002 


\CroV 






AG220 = AG020 + AG200 


[Cro]^ 


Zx 




AG202 = AG200 + AG002 


[Cro]^ 






AG120 = AGioo + AG020 


[cl] [Cro] 


Zx 




AG210 = AG200 + AGoio 


[cl] [Cro] 


Zx 




AG102 = AGioo + AG002 


[cl] [Cro] 






AG201 = AG200 + AGooi 


[cl] [Cro] 






AG012 = AGoio + AG002 


[cl] [Cro] 






AG021 = AG020 + AGooi 


[cl] [Cro] 






AGiii = AGioo + AGoio + AGooi + AGcoop 


[cl]' 






AG222 = AG200 + AG020 + AG002 


[Cro]' 






AG112 = AGioo + AGoio + AGcoop + AG002 


[cI]'^[Cro] 






AG121 = AGioo + AG020 + AGooi 


[ciyiCro] 






AG211 = AG200 + AGooi + AGoio + AGcoop 


[clYlCro] 






AG221 = AG200 + AG020 + AGooi 


[cI][Cro]^ 






AG212 = AG200 + AGoio + AG002 


[cI][Cro]'^ 






AG122 = AGioo + AG020 + AG002 


[cI][Cro]'^ 





Table 2.1: Values of variables referenced by the partition function. 



^mpi is phages per cell. 

values for these states are also actively measured by the experimental community 
and could be included in the partition function insofar as RNAp does indeed bind 
in the Or, so that in any given interval of time there is indeed the probability that 
the Or is in a state not considered in the above partition function. The effect is 
said that eliminating the RNAp is constant, though there is a continuation in the 
literature that prefers to keep RNA polymerase explicitly in the partition function 
and no mention to what the value of that constant. The situation is a little subtler 
then that and the value of this effect, as shown here, relates to the RNAp binding 
energies to the Pr, and Prm as well as absolute concentrations. The benefit is 
still a simpler partition function and since RNAp is fixed for the simulation, the 
additional monomials do not alter the algebraic structure of the equations but only 
append to the weights of the already present monomials, but this effect can still 
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induce bifurcations and one should be careful. 

It is then necessary to write the dimer concentrations in terms of the monomer 
number. Here we use the expression implicit to the results of Zhu et al. [2004], 
Wang et al. [2010], where we explicitly show how the number pnc, labeled as a 
conversion factor from number to concentration, enters into their equation. 

We note that the above equations appear with the mistake of a positive sign, in 
Wang et al. [2010] and the arxiv version of Zhu et al. [2004], and use the expression 
from Reinitz and Vaisnys [1990] with the insertion of the pnc, a converting factor 
from protein number to concentration, and implicit assumptions about RNAp 
(see the note on dimerization). The bacterial volume used in Zhu et al. [2004] is 
stated to be .7 * 10~^^L and pnc is in inverse proportion to the volume. The value 
used for pnc is 1.5 * 10~^^. A calculation reveals and is corroborated by a private 
communication that in actuality pnc = {volume * Na)~^ = 2.37* 10~^ Moles/ L 
and the listed value is a typo. Use of that number reveals unrealistic protein 
concentrations and it is this authors belief that the number pnc includes implicit 
assumptions about the RNAp concentrations and affinities. Provided RNAp is 
removed from the partition function, a pnc value of 1.5 * 10~^^ is indeed more 
consistent with the numbers obtained from experiment, and was likely tuned ad 
hoc to the present value for that reason. No reference is made about the value of 
the number or conversion in Wang et al. [2010], although the calculations present 
suggest the use of the conversion from Ao with the equation from Reinitz despite 
the errors in their written equation. Many calculations should be performed twice, 
once with the values similar to Zhu et al. [2004] and again with RNAp in the 
partition function and pnc = 2.37 * 10~^ which is postponed for future work. As 
noted in the introduction this equation should be modified in future work anyway 
as a result of the higher order polymerization processes that occur. 

2.2 A Partition Function; RNAp or no RNAp? 

The partition function in Shea and Ackers [1985] reads, with [RNAp] the concen- 
tration and assigning the index ^3' to denote the RNAp occupancy. 
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Z = 1 + ([RNAp](e-^^o3 + g-^Ggo) ^ [RNApfe-^^^^) + 
[c/] + g-^^oio ^ g-AGioo _^ _^ g-AGios _^ g-AGi3o)^_^ 

[Cro] + g-AGo20 ^ g-AG200 ^ [RNAp] (e"^^32 + g-AG203 ^ g-AG230))^ 

[c7]2(g-^Goii ^g-AGioi ^g-AG„o ^ [R]VAp]e-'^^ii3) + 

[Cro]2(e-^^022 + g-AG202 + g-AG220 + [R]VAp]e-^^223 ) + 

[c/] [Cro] (g-^'^on _^ g-AGo2i _^ g-AGio2 ^ g-AG2oi _^ g-AG2io _^ g-AGi2o_^ 

+ [RNAp] +e-^G2i3))^ 
[c7]3(e-^Gni) + 
[Cro]3 (e-^G222) + 

[c7]2[Cro](e-^^ii2 + g-AGm + g-AGzn)^ 
[C/] [Cro]2(e-^Gi22 + g-AG2i2 + g-AG22i) 

(2.7) 

It is clear that the effect of [RNAp] is not constant; [RNAp] does not factor out 
of the partition function. Three experimental numbers, the RNAp concentration, 
the RNAp binding to the cl premoter, and an RNAp binding to the Cro premoter 
are entangled into this partition function. Removing RNA polymerase from the 
partition function also has the effect of disturbing the time scales, in that equation 
2.1 only accounts for the fraction of a dt in which no RNA polymerase is bound. 

We note that this form of the partition function would require corrections for 
low numbers in that it uses the approximations ^ n^{n — l),n^ ^ n{n — l){n — 
2). It is likely that some of the complexity of these systems near the axis is residual 
of such an approximation. This becomes especially more important if multiple 
phages have invaded the cell or the Oi is included in the partition function. Implicit 
to the figures in Zhu is the T^, Tmif and Trr values reported are 1/3 of the values 
used in the calculation. This is due to the inclusion of 3 phages in the cell as in 
Aurell and Sneppen [2002]. The calculation left out presumes that = Z^, and 
the states that promote a ci, for example, are enhanced by the number of phages 

present . Also, it is often presumed that for the composite system of the left 

and right operator regions Z^^ZqiZor- For a partition function to be the product 
of the subsystems, the subsystems must not interact. This is not the case here, and 
especially for the low numbers involved. One need only plot, for example, next 
to X * (x — 1) * (x — 2) * (x — 3) * (x — 4) * (x — 5) * (x — 6) * (x — 7) * (x — 8) to see 
that the difference is considerable. This error is further increased with monomers 
present as well, since this too reduces the number of dimers and increases the 
fluctuation. 

We use AGci = —11.1, AGqw = —7.0, and AGcoop = —6.9. The states and their 
values are collected in Table 2.1 
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The values for the multiply liganded states are found by summing up the 
values of the constituent singly liganded states. Cooperativity gives an additional 
sum, AGcoop = —6.9, for the states with cl and cl on adjacent sites. These are 
the G*ii, Gi*i configurations where * can denote 0,1, or 2. RT = .6V7kcal/mol is 
the Boltzmann factor for the 310 Kelvin laboratory conditions which these were 
accumulated. When we speak of the mutant A/yj^ we mean that the affinity for the 
ith site is used in place of Or^, jth in place of Or^, and kth in place of Or^, Some 
of these are likely biologically unstable but it is inevitable that different strings for 
the operator configurations may have similar AG and thus similar dynamics. 



3 Dynamics Near Equilibrium 



The key to the treasure is the treasure. 
John Barth, Chimera 

3.1 Gene Expression and General Dynamical Systems 

Although this is a simplified model, the need for general dynamical systems is 
already evident. 

The number of stationary points is complicated by dimerization but in the basic 
case this number is found to be in proportion to the {number ofbindingsites + 1) x 
(numberofproteinspecies). The positivity of the partition function and concentra- 
tion ensures that some of these are physically realizable. The dimerization used 
here doubles the number of stationary points. Thus it is reasonable to start from 
very general mathematical considerations Liapounoff [1893], Perko [2001] and well 
over a century of elegant mathematics that would be required to understand the 
mathematics of gene expression. One is compelled to introduce a modified model 
that is more general and symmetric in the form in which creation and annihilation 
processes enter. For example, consider the n equations to replace for the discussion 
equation 1.3. 

Vxi^x^ (xi, Xn) H2O, RNAp, Rihosomes, ...) — y^\^^\ Xn) H2O, RNAp, Ribosomes, ...) 

— i i 

^ Z(xi, Xn) H2O, RNAp, Rihosomes, ...) 

(3.1) 

Here one can begin with a partition function for the n proteins and let quantities 
like H2O, RNAp, all the things that contribute non-negligible partial pressures in 
the cell, enter as control parameters. Say, Z{xi,...,Xn)H20,RNAp,Ribosomes,...). 
Again Vx- is the rate at which transcription in the Zx- states occur. Introduced is an 
expression for annihilation, v^tZ^t, to account for the specific events that destroy 

a protein. For example, in the A phage situation, the cl protein is destroyed by 
another protein, RecA, which destroys a cl monomer by cutting it in half. In such 
a situation one would include a term Z^jt oc [cl] [RecA] governed by a rate, v^jt of 

this destruction process. 

The motivation is that the stationary points would appear at the points Xi, 
where VxiZx^ (xi, x„; H2O, RNAp, Rihosomes, ...) — v^t Z^t (xi, x^, H2O, RNAp, Rihosomes, ...) = 
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for all i. This is a set of polynomial equations in the protein variables X/. This allows one 
to either fit to known systems, or construct others. Also since, Zx^y^z is positive, and any 
term appearing in the numerator appears with a positive sign in the denominator, many 
qualitative features persist. Say x = /(x) and y = f{y)/g{y)' One can see that if x = 0, 
then y = 0, since g(y) is non-negative that the stationary points in the first system are the 
same as in the latter. Next, in the neighborhood of a stationary point the dynamics are 
qualitatively the same as the linear system. 



dx 



^stationary) / (3-2) 

^stationary 



and 

9/ 



|^(y) - i/(y) 

^ - 



g{y) 



2 



(y y stationary) ' (3-3) 



y y stationary 

Since at a stationary point f{y) = 0, evidently. 



a/ 



8iy) 



(y y stationary) ' (3-4) 

y—y stationary 



The two systems have similar quantitative behavior then too, with the dynamics simply 
scaled: 

t^t^giy) . (3.5) 

y—ystationary 

For an overview of such methods in dynamical systems many texts are available, Perko 
[2001]. 

This interlude is to motivate some of the extent to which these kinetics serve as 
analog circuits. Restricting the numerators to polynomials, //(xi,..., x^)/ in which fj is 
non-negative when Xj = 0, that is to say we require proteins to exist before we can destroy 
them. Then we can always program a set of binding-expression interactions to simulate such 
dynamics. 

For example, say one would like to construct a circuit with dynamics similar to 
x = x + x^ — x*z, y = y + y^ — x*y, 2; = z + z^— y*z. Weighting by Zx,y,z = 1 + x + x^ + 
y + y'^ + z + z^ + x^y + x^z + y^z does not change the locations of the stationary points. 
That is to say, one can accomplish this circuit with three proteins, with two binding sites 
each, in which the occupation of either would enable their own transcription; also the z 
protein destroys x, y protein destroys z, and x destroys y. 

To be able to fully fit the experiments of St. Pierre [2009] the phase space would 
partition into more than just two states, but four states, in addition to a splitting of a cell 
population into lysis and lysogeny, some cells (refractory) were immune to the phage and 
a small fraction of the lysogens had daughter cells that lyse. The number and properties 
of lysogens that had daughters that underwent lysis (5/100) could be useful in additional 
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studies of robustness ^. We continue the analysis of Zhu et al. [2004], Wang et al. [2010]. 
Including the Ol, RNAp, tetramers, etc. could contribute this additional complexity, codes 
are provided Borggren [2011] if one was inclined towards such study. 

At different times in the cell cycle different genes are turned on or off. Often these 
genes are transcripts for proteins that are necessary for functions in that particular interval. 
The concern here is with the lysogeny decision and in this interval we have only the Or, cI 
gene, and Cro gene on. 

Letting x, y denote the cl and Cro monomer numbers respectively we write for the 
dynamics 

X = VxZx{x,y)/Z{x,y) - x/x^, (3.6) 

y = VyZy{x,y)IZ{x,y)-ylXy. (3.7) 

The creation term in the dynamics of the associated protein populations within a cell 
are governed by the rates at which RNA polymerase binds and translates the DNA into 
mRNA, together with the rate at which the ribosome translates the mRNA into the chain 
of amino acids that comprise the protein. In this way a cell seems to possess three 
different molecules encrypting the same message. Categorically speaking, the fundamental 
biological information processing is governed by the physics of RNAp : DNA RNA, 
Ribosome : RNA^ amino acid. The net effective rate of these processes enters directly 
into the equations of motion through the Vx,y variables. 

One can say much about these questions without specifying numbers of binding sites 
or which promotional states enter into Z/(x,y) if we grant a few specifications. Firstly, 
we assert the positivity of and i/^. This is also largely natural selections choice, though 
alternatives are imaginable. For example, RNAp can indeed unscribe a transcript, the 
backward reaction is permissible, and it is said that for any 11 forward reactions, there 
were 10 steps back Bennett [1982]. The partition functions are positive by virtue of being 
well-behaved non-negative partition functions or a subset thereof. To avoid singularities 
and to insure Xi > we include at least one state that does not depend on x and y in 
each of Z, Zx, and Zy. In the case of A this is provided by the open DNA configuration. 
This is the reference state AGqoo = 0- In the original formulation. Shea and Ackers [1985], 
security against singularities was also provided by the explicit inclusion of [RNAp] and 
[RNAp]^ states. For two-dimensional systems, for example two genes expressing two 
proteins, this is sufficient to prove the existence of any non-negative number of stationary 
points provided a sufficient supply of binding sites. 

Stationary states occur when the derivatives vanish: 

xZ{x,y) - VxrxZx{x,y) = 0, (3.8) 

yZ(x,y) - VyryZy{x,y) 0, (3.9) 

^ An interesting explanation without requiring more basins of attraction would be, if, for example, 
as a cell divides the daughter cells have each taken half the cl dimers from the parent lysogen, and 
in so doing switch to the lytic basin of attraction, causing the cells to lyse. (It is imaginable that 
whereas x, y may sit in one basin of attraction x/2, y/2 may reside in the other) 
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This is to illustrate some of the most important aspects of the equation and how these 
states interplay to create the desired states. 

When X, y are zero, x,y are positive since Vx^y and Zx{x,y) / Z{x,y) are. As x, y go 
to inf the —x/Zx and -y/ Ty terms dominate and since Tx^y are positive, the dynamics are 
confined to the positive quadrant. It follows similarly for y as well as x that there exists 
a contour for which x = on their respective axis since < Zx{x,y) /Z{x,y) < 1; thus x 
switches eventually to negative and stays negative and this crossing requires at most VxTx 
proteins. 

This provides an experimental test, or at least an avenue into parameterizing the 
models. For example, if one were to measure only Vx^y, Tx^y, and just the expectation value 
of Zx,y{x,y) / Z[x,y), one already has means to predict the average number of proteins x 
and y. Since the expectation value of Zx^y{x,y) / Z{x,y) is just the fraction of time RNAp is 
promoting, it is a much easier process than discovering all possible states, measuring their 
relative affinities and then simulating. 

It is the fitting of empirical models to experimental data that comprises the bulk of the 
quantitative features of biophysics. 

With cl at 236 amino acids, requiring at minimum 708 base pairs and Cro, 66 amino 
acids for 198 pairs, it is not surprising that Cro can be transcribed faster. The decay time 
in seconds is Td — 2943 and Xcw = 5194. 

In this model the states that contribute to the numerator for x are those in which Oro, 
is unoccupied and thus open for RNAp binding. Of these states those with a cl dimer 
bound at 0^2 ^Lre believed to further enhance the rates. Thus, 



cl _1 

'7 \ 



dt 



mpz * Trm * Eci{[cl\e ^ + [clYe ^ + [cJ][Cro]e ^) 

+ Trmu * Ea(l + [cl\e ^ + [Cro\{e ^ +e ^) + [cl\[Cro\e ^ + [CroYe ^) 
. cl 



(3.10) 



Cro _ 1 . 
lit ~Z^ 



mpi * Trr * Ecro * (1 + [cl\e + [Cro]e rt (3.11) 
Cro 



) 



Here Trm/Trmu/ airid Ed are 0.115 cl monomers/sec, 0.01045 cl monomers/sec and 
1 copy per mRNA synthesized. Next, RNAp is assumed to require both Ori and Or2 
free to bind to the Cro promotion site. This occurrence is repressed in comparison with 
the states that promote cl but promotions proceed with a higher burst. Trr.Ecto are 0.30 
Cro monomers/sec and 20 copies per mRNA synthesized, respectively, mpi denotes the 
number of phage Or sites per infected cell. The mpi factor curiously seems to only have 
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been multiplied by the premotional states with cl dimer at Or2 arid has been given the 
value 3. Here, Z is the partition function from 2.1, while [cl] and [Cro] are the assumed 
relations in 2.3 and 2.4 respectively. Again, these are the assumptions used to correspond 
to Zhu et al. [2004], Wang et al. [2010] and to references therein. We will posit different 
equations of motion, including the correction to the pnc confusion, and the effect of RNA 
polymerase in variations on a theme. Experimental collaboration, references, or data to 
analyze that can pin down any to all parameters are welcome, appreciated, and daresay 
necessary. ^ 

The situation is more complicated in that cl promotion occurs at 2 different rates, the 
states with a cl dimer on Or2 enhance the actual rate which RNAp transcribes, as well as 
binding more favorably, Ptashne [2011], and how this assistance translates physically into 
the enhanced rate is an interesting open question. Here it is expressed by using the two 
different transcription rates, Trm^^rmu- 

To approximate the effects of specific mutations we consider mutations, Xip by which 
it is meant that new cells have been synthesized with Or\ Ori ,Or2 Orj and 
Or3 Or]^, the chemical mutations correspond to changes in the strings shown in Fig. 1.4 
which lead to changes in the binding affinities in Table 2.1. 

3.2 Stationary Points 

The stationary points of the system are found empirically and tabulated herein. 

To inquire into the nature of these stabilities it is fruitful to consider the linearized 
equations of motion at a stable point which approximates the flow in the neighborhood. 
Stable, unstable, and saddle points have negative, positive, and mixed eigenvalues for 
their respective linearized equations. An imaginary part suggests the dynamics around 
the point have non vanishing circulation. 

Although we will maintain a global perspective in the Diffusion section to follow, 
in correspondence with the main utility of the methods in Wang et al. [2008] or Wang 
[2011] where a local analysis is not necessary before proceeding to find global features, 
dynamics near stationary points is mathematically very well understood Perko [2001] 
and even with stochastic features Kwon et al. [2005], and any such analysis would begin 
with an enumeration of the fixed points provided here. It is also noticed that the peaks 
of the distributions to be found are displaced from the fixed points of the underlying 
dynamical system. This displacement is purely the result of the inclusion of noise and 
has been observed in theory and numerical studies before Kwon and Ao [2011]. An 
experimental situation where the classical dynamics can be fixed and the noise part can 
be scanned independently of these dynamics could greatly facilitate the observations and 
interpretations of the non-equilibrium steady states and the breaking of detailed balance 
suggested in Wang et al. [2008], Wang [2011], Qian [1998], Ao [2004] by directly measuring 
this displacement. ^ 

^nborggren@gmail.com 

^ Doing this in a minimally biased fashion could be difficult. For example, scanning noise by 
simply letting the number of particles decrease can easily put one in a domain where the classical 
kinetics are no longer a good approximation Samoilov and Arkin [2006]. In such an instance 
the measured displacement, in addition to the contribution from noise, would also include a 
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For these systems the relevant matrix components are, writing Xi,...,x„ as x, 

Vf{Xi,...,Xn)ij = ^^^^-^^^^^[Z(xi,...,X„)^Z;,.(Xi,...,X„) 

- ZxjiXi, Xn) g^2;(xi, x„)] (3.12) 

/ 

where 5ij is one when i = ; and is otherwise zero while the partition functions are evaluated 
at the stationary point under query. 

For the two by two case the eigensystem is 

Oil = - lad + ^hc + (P- + a + d), (3.13) 

(3.14) 



/ VT^^ - 2ad + Ahc + d^)-a + d , 
^1 = ( Yc 



0C2 = -(Vci^ -2ad + 4bc + d^ + a + d, (3.15) 



/ VT^^ -lad + 4:bc + d^)+a-d , , 

^2 = Yc 



(3.16) 



with 



c=|[zAz,-z4zl, (3.19) 

In two dimensions the quantity {ad — be) is the determinant. The fact that this term 
need not vanish is why the notion of energy landscape must be generalized into the flux 
framework to include fields that are not the gradient of a potential, see Wang et al. [2008]. 

Understanding these equations can tell us a great deal about the type of forces at work 
since near a particular point, xf = (xi/, x^/) through an expansion Perko [2001], we 
know that the system is qualitatively similar to the linearized equations 

x^ f{xf)+Vf{xf)x. (3.21) 



contribution from reality deviating from the kinetic assumptions of the model. 
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If Oil <0, 0C2 <0 then the point serves as a probability sink, if oci > 0, 0C2> 0, then it is 
a probability source. Mixed eigenvalues are transition states. For example, one can verify 
that any of the following conditions lend to stability: 
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0, 


d < 







a = 


0, 
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1^-bc. 



Of the 27 mutations considered, 15 had underlying dynamical systems with 3 stationary 
points, 11 had one stationary point, and one had two stationary points. The locations are 
tabulated in the tables A.l, A.2, A.3 and A.4. In particular we note that the permutation 
group S3 acting on Ori,Or2, and Or2>, only six mutations, is sufficient to provide mutants 
with a single lytic peak A231, A321 a single lysogenic peak A132 while Ai23,A3i2, and A213 
retain the switch properties. 

The 27 phase portraits of this model can be seen in Tables A.l, A.2, A.3. Dynamics are 
also combined with the results of the diffusion section to demonstrate the importance of 
statistical, trajectory based approaches. Those results are in the tables A.25, A.26, A.27 and 
demonstrate that the different limit cycles correspond to regions of different entropy or 
energy density. 

This model, however, does not probe what the downstream reactions do to the phase 
space distribution given to it; it is quite likely that in the case of, for example, A231 lysogens 
will still attempt to form but with a loss of stability. It does, however probe transient 
phenomena and the initial impulse of the system is, even in the monostable case, towards 
two states. 



dZx ^\cl\ , AC 



dx 



dx 



010 _\_ p-^Gioo\_^ 



2*[cJ]*M,-AGno + 

d[cl] 



dx 



dx 



9y 5y ' 



.d[Cro] 



-AG220 



+ 



2[Cro]^^e 
dy 

[cJ]^[g£l(e-^G2io+,-AGi2o) 



(3.22) 
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The robustness of the lysogenic cycle would be unsurprising if the system were wholly 
deterministic. Dynamical systems in general can partition into multiple basins of attraction 
and thus so too can those arising from chemical kinetics. In two dimensions, like that 
presumed here, it can often be the case that there exists some contour that partitions the 
phase space Perko [2001]. This contour is unstable in the perpendicular direction and 
deterministic trajectories beginning on this contour take the steepest descent to a saddle 
point if the system is bounded. Any crossing from one region to the other would require 
an intersection with this line and is forbidden with the rules for ordinary differential 
equations. Thus, hopping would require a stochastic process to change this fact for this 
model. The puzzle of this stability, which we consider in the next section, is why there is 
not enough noise. 



4 



Diffusion; Evolution of Ensembles 



Lo otro no existe: tal es lafe mcional, la incurable creencia de la razon humana. 
ldentidad=realidad, como si, a fin de cuentas, todo hubiera de ser, absoluta y necesaria- 
mente, uno y lo mismo. Pero lo otro no se deja eliminar; subsiste, persiste; es el hueso 
duro de roer en que la razon se deja los dientes. 

Antonio Machado, ]uan de Mairena, quoted in El laberinto de la soledad, 
Octavio Paz 

When Galileo first assisted his gaze with the telescope and directed it towards Saturn, 
he discovered its ring, and Titan, one of its moons. The point-like orb of the ancients 
revealed more structure still when Cassini subsequently fixed his gaze upon this ring, 
revealing a large division; the ring was not one, but two. Further query revealed yet more 
structure, courtesy of an in situ spacecraft, also named Cassini, and one is now disinclined 
to the naming of rings, which appear with the abundance of grooves in a phonograph. 

Indeed the rings of Saturn are not rings at all, but a collection of granules of various 
sizes caught in circular potential wells; the only rings are in the energy landscapes, our 
abstractions, and only appear in the data in that there was a sufficient supply of granules 
to sample the space so that the typical distance between stones was shorter than the 
resolution of our instruments. 

With that said let us revisit the diffusion equation 



For, example if one were to construct a Hamiltonian of the Saturn system and fill 
a vector, F^, with each of the equations of motion, then u would denote the density of 
granules, and in the case of Djj being zero, the equation would describe the evolution of 
that initial configuration as it samples the Saturn force field; depending on the details 
some fraction fall into Saturn's core, some might fall into Titan's surface, others escape, 
and some get trapped into the resonances of the system, the rings. Letting the D > 
would allow our granules to interact with each other, although they are small enough so 
as not to perturb Ff, they are of comparable sizes to each other and collisions lend them to 
diffuse through space. 

In the case of A phage, the numbers of proteins involved are small, and intrinsic noise, 
shot noise, is indeed a real source of diffusion for the system; the system kinetics do not 
behave the same with every trial, but spread throughout the phase space. However, Eq. 
(4.1), can equally describe such a situation. For references and the general state of affairs 
of Fokker-Planck equations used as a theory of everything, query Friston and Ao [2012]. 



dfU = 



di{DijdjU — Fju). 



(4.1) 
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4.1 Master Equations, Langevin Dynamics, or Diffusion Equation 

Many approaches are used to capture the chemical kinetics or gene expression in a 
noisy environment, all of which are phenomenological and have their advantages and 
disadvantages. The proper approach can vary from system to system and depend on 
resources available. See for example Wang et al. [2008], which this thesis follows and 
extends with the finite elements, and also the potential approaches in Ao [2004], Zhu et al. 
[2004]; for field-theoretic approaches there is Wang and Wolynes [1996], Wang et al. [2010], 
Feng and Wang [2011], and for Langevin dynamics and master equations see Wang et al. 
[2006], Ge and Qian [2010], Qian et al. [2002]. 

4.1.1 An Integral Approach for Diffusion Equations and Passage Times 

The finite element method is emerging as a universal method for the solution of differential 
equations of arbitrary complexity from all walks of science, see Logg et al. [2012]. The 
following derivation is typical of the procedures implemented throughout the FEniCS book 
(Logg et al. [2012]). In order to use a finite element method to solve partial differential 
equations it is necessary to first translate the equations into an integral equation. This 
integral equation can be seen as setting a bilinear functional equal to a linear functional 
and solving the corresponding system by standard linear algebraic methods. It will 
be described in detail for the Fokker-Planck equation and corresponding passage time 
equations. 

4.1.1.1 Diffusion Equations 
First we seek to solve: 



For numerics the function space is specified by defining a mesh over O and the type 
of element associated with each point in the space. These specifications alter the analysis 
mainly in how long it takes for a desired level of accuracy to be reached. We continue 
thinking of the function space as the space of twice differentiable functions over O. The 
functional form is found by multiplying a test function v and integrating over the mesh 



dtU = di{DijdjU — Fiu). 



(4.2) 




(4.3) 



Continuing with integration by parts yields 



v{DijdjU — Fiu)nids — / div{DijdjU — Fiu)dx. 



(4.4) 



Here D/y denotes the spatial dependance of the Diffusion, F/ , is the drift force, n/, a normal 
vector to the mesh. The mesh is denoted O with boundary 90. 

We approximate to find an expression for Uf+st given the solution at Ut 
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/ v(ut^st — Ut)dx ^St[ / v(DijdjU — Fiu)nids 
Jn Jdct ^ ^ 

— J div{DijdjU — Fiu)dx]. 



(4.5) 



We take the average of the expression over time St hy asserting u ^ Mi+^±lii Continuing 
gives 

St f 

y [ J^^v{DijdjUt+st - FiUt+st)nids 
- / div{DijdjUt+st - FiUt+st)dx 

y (4.6) 

+ / v{DijdjUt — FiUt)nids 
div(DijdjUt — FiUt)dx]. 



Collecting terms with Uf^st ori one side of equality and Ut on the other gives the variational 
forms a and L that we will use 

f St 

ci = J^(^^t+^t + —div[DijdjUt+st - FiUt+st])dx, (4.7) 

L = J^i^ut - ^div[DijdjUt - FiUt])dx. (4.8) 

a{v,u) = L{v) then defines a linear system. To see this note that if v, v' satisfy a then so 
does oc{v + v') and similarly for L for all scalars oc. Since it also holds that a is linear in 
UtJ^stf a is bilinear. The Ut in L is specified by propagating uq and thus explicitly appears 
as a known function on the right side of the equations. Thus, finding the u that solves 
a{v,u) = L{v) for all v advances our solution a timestep. Eq. (5.1) shows the simplicity of 
the lines of code to perform these calculations. Indeed the computation is more general 
and concise then the notation here. 

The passage time equations are much simpler. The integral equation corresponding to 
Fiditi{x) — Ddiditi{x) = — 1 follows by simply multiplying through by the test function 
and integrating 

/ {t{x){Fiditi{x))-Ddit{x)diti{x))dx = - [ t{x)dx, (4.9) 

[ (t(x)(Fiditn{x))-Ddit{x)ditn{x))dx = -n [ t{x)Tn-idx. (4.10) 
Jn Jn 

One can then find the passage times recursively. The induction step is bilinear provided 
one has first solved for T^_i(x) before proceeding to find Tn{x). 

To perform the simulations in this work we need still to specify a mesh and type of 
finite element. We choose for O the domain in protein number space < cJ < 1000 
and < Cro < 1000. We discretize this into at least a 50 by 50 grid and use a quadratic 
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Lagrange element. The quadratic Lagrange elements span the surface by connecting the 
adjacent mesh points with a quadratic function (see Fig. 4.1). The function space is the set 
of all such functions. One benefit of a finite element method over a finite difference is the 
functions in the function space are continuous over O and we can evaluate any function, 
and in particular our solutions, anywhere in the domain. We see this smoothness manifest 
in the plots which evaluate the functions along the trajectories. 

Figure 4.1: Finite 
elements construct 
a surface by com- 
bining simplices 
together. The 
quadratic elements 
used fit the func- 
tion inbetween the 
mesh points with a 
quadratic formula. 




4.1.2 On Convergence 

Technically, the asymptotics of the Fokker-Planck equation take infinitely long to converge. 
Such a program, see Fig. 4.2, can still be written, albeit not technically admissible into 
von Neumann architectures. Taking the a and L as bilinear and linear forms, respectively, 
defined anywhere, we have an infinite trajectory pursuing eternally its asymptotic states. A 
break can be entered by establishing a criterion for convergence. The distributions present 
in the Appendix were found using 53, 7.5 second time steps and well approximate the 
steady state solutions. 

This loop was carried out 53 times before I choose to break it in the states labeled 
steady in this work, with k = 7.5 'seconds' the time step. 

4.1.3 Choice of Gauge 

It is noticed that equation (4.1) presupposes a choice of gauge and the curl of any vector 
field can be added to the flux landscape, DjjdjP — FjP + e^^^djA]^, and still satisfy the basic 
equation. We choose a gauge such that e^i^djA]^ = 0. We notice that what has been called 
the entropy production in Wang et al. [2008] is a gauge-dependent quantity and how it 
translates into entropy or heat in absolute units is open. We note that the different kinetic 
paths will have different thermodynamic signatures, and understanding the connection 
can enable direct measurement of some of the quantities that have been proposed. 

A recent paper has elaborated on the nature of the relationship between related Fokker- 
Planck equations and gauge theories, Feng and Wang [2011], this is likely to prove a rather 
pregnant approach to chemical kinetics. 
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Figure 4.2: This 
would be an exam- 
ple of a legitimate 
python code that 
runs indefinitely. 
Computationally, 
the asymptotics 
of the equations 
require a function 
that does not 
terminate. 



Python code 



def TimeStep(u0, ul) : 
A = assemble(a) 

while 'no your solutions have not converged' 
b = assemble(L) 
bc.apply(A,b) 
solve(A, ul.vectorO, b) 
u0.interpolate(ul) 
yield ul 

while InfiniteTraj : Infinite! raj .Next() 



Figure 4.3: Calcu- 
lating expectation 
values is straight- 
forward, even on a 
command line. 



Python code 



»>exprs = [Expression! "x[0] "), Expression ( "x[l] ") , 

. . Expression ( "X [0] *x [0] ") ,Expression("x[l]*x[l] ")] 
»> ans=[ 

assemble(i*u*dx,mesh=mesh)/assemble(u*dx,mesh=mesh) 
for i in exprs] 
»> ans 

[93 . 076476131804114, 169 . 77672009751805 , 

13624 . 244707344435 , 36974 . 455745161322 ] 
»> sqrt (abs (pow(ans [0] ,2)-ans[2])) 
70.434468112069482 
»> sqrt (abs (pow (ans [1] ,2) -ans [3] ) ) 
90.279128585129428 



4.1.4 AxA/ 

After having computed the solutions to the equations we can calculate some moments of 
the distribution. So far the diffusion has been picked to correspond to the numbers believed 
to quantify robustness, see Zhu et al. [2004], Wang et al. [2010], however the boundary 
condition can be tuned to make this number vary over many orders of magnitude and it 
may be more prudent to use moments of the distribution to parameterize D. For example 
when cl proteins are reported in numbers 150-250 in a cell this can suggest an expectation 
value of, say, 200 for x and a spread of 50 for the second moment. These numbers are 
more easily measureable then the robustness, which has been the subject of amendment 
over time. 

The figure 4.3 illustrates how FEniCS software can compute expectation values from the 
command line, and the type of calculations that can be performed once one has solutions 
in hand. Let us compute Ax and Ay with the distributions from A.IO and recompute the 
related expectation values. 

It is important to realize that it is not the mass of an electron that biologists measure. 
The errors do not necessarily go away with better measurements. The numbers are 
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supposed to fluctuate and this can actually be interpreted to mean that the higher moments 
of the distribution are relevant. A spread of 70 in the number of cl proteins and 90 Cro 
proteins is predicted by this model. 

4.2 Passage Times 

The equations again for the first passage time is 

/ {t{x){Fiditi{x))-Ddit{x)diti{x))dx = - [ t{x)dx, (4.11) 
Jet Ja 

where we force a boundary condition of zero, in the present case a circle of radius ten 
is placed around the critical point, where we want to calculate the passage time to and 
following we tabulate the results for the mutants we can construct from the binding values 
available. First, we explore the properties of the 15 of 27 mutants that maintain the switch 
properties in this model. We notice however that the illusion of bistability is manifest 
in the transient phenomena even for the monostable states. That is to say, if the switch 
behavior is a non-equilibrium phenomenon, monostable mutations may be relevant. 

4.3 Noise-induced stability 

The dynamics picture of the previous chapter is quite idealized. Even if all the parameters 
were perfectly known there is still the issue that we are discussing mesoscopic regimes 
and our large N intuitions and classical assumptions are certain to break down in this 
limit. The constant rate assumptions would need to be replaced with discrete events and 
the manner in which a unit of time splits into subunits becomes important, as well as 
the annihilation processes. The changes to the partition function in this model are again 
relevant. 

It is found here that the noise added displaces to a large extent the locations of the 
peak from the corresponding stable points of the underlying dynamical system. Viewing 
the noise as diffusion through concentration space and developing effective field theories 
and path integrals has provided interesting and useful constructions, see Wang et al. [2010]. 
The fact that there even exist corrections to be made is interesting. This is the type of 
functionality that Schrodinger deliberates on being prohibited for an organism, in his 
discussion of the \/\n) law. 

... an organism must have a comparatively gross structure in order to enjoy 
the benefit of fairly accurate laws, both for its internal life and for its interplay 
with the external world. For otherwise the number of co-operating particles 
would be too small, the Taw' too inaccurate. The particularly exigent demand 
is the square root. For though a million is a reasonably large number, an 
accuracy of just 1 in 1,000 is not overwhelmingly good, if a thing claims the 
dignity of being a 'Law of Nature.' 

Erwin Schrodinger 'What is Life' 

For A phage this is just a restatement of the stability puzzle: how then can the numbers 
of proteins be within the range 0-1000, a far worse scenario then Schrodinger deemed 
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Figure 4.4: The 
first passage time 
with an absorbing 
boundary condi- 
tion within 10 
units of the origin. 
The origin is ex- 
tremely unstable; 
The steepness 
of the funnel of 
the first passage 
time landscape 
illustrates the 
overwhelming 
propensity for 
the system to be 
driven towards 
some non-trivial 
state. The num- 
bers, ^ 10^, appear 
in seconds. 



Python code 



Assembling matrix over cells [======================== 

==============] 100.0% 

Applying boundary conditions to linear system. 
Solving linear system of size 1442401 x 1442401 

(PETSc LU solver, umf pack) . 
*** Warning: Using LU solver, ignoring preconditioner 

"default" . 

evaluate at wells, [0.14175299999999999, 

462 . 57299999999998] 112474076 . 373 
evaluate at wells, [644.05200000000002, 

. 0061165799999999999] 112473547 . 501 
evaluate at wells, [29.979800000000001, 

178 . 22999999999999] 112469292 . 141 



om 




possible? Why does a deterministic dynamical analysis in the small concentration regime 
work at all? Do we retain the observed robustness when we consider noise? Is the noise 
really vTn) noise? What paths are permissible and what are their relative probabilities? 
General questions and methods of this variety are the subject of intense debate and muse 
to much creativity, see Wang and Wolynes [1996], Wang et al. [2010], Stock et al. [2008], 
Ghosh et al. [2006], Qian et al. [2002], and have been applied to the A system, see Zhu 
et al. [2004], Wang et al. [2010]. This debate has led to legitimate concerns and elegant 
measurements to seek explanations from other variables, like volume, where the variation 
across an ensemble are less severe or can be experimentally minimized, see St. Pierre 
[2009]. 

It is written that DdCro = const * TciCro/^ciCro then A = D~^, Zhu et al. [2004], Wang 
et al. [2010]. However, their expression for A is the same as their expression for D. The 
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calculations require 

/ 17.85*^ \ 

644.052 corresponds to the value at the lysogen zero and Ncro = 462.573 the 
Zci = 2943s and Zcro = 5194s giving 

°»' " ( 'o' 2.2271 ) ■ 

The numerics proceed with this diffusion tensor. If the analogy is to extend from diffusion 
in ordinary space to diffusion in concentration space, then the mean squared displacement 
is a displacement in concentration space and the units of this D should be dNcidNcw/ds. 
Let us attempt to motivate the form of Eq. (4.12) and the scale of Eq. (4.17). If 

X f ( x^ 

X ^ x±^/x then x ^ x± — ^ ^ fix) ± (4.14) 

2^/x ^ ^ ly/x 

For large x, we have seen that the term —x/r dominates and recover fluctuations of order 
^/x. Since D goes as the square of these fluctuations we get the proportionality to x 
as appears in 4.12. The situation is more subtle when x is small. For Shea and Acker 
type models the Dxx term, say, would be dominated by terms of the form Z^/ xZ^ which 
can diverge. We expect corrections from this end to increase the average diffusion and 
therefore facilitate the crossing of a barrier, thereby reducing the robustness of the switch 
more than Eq. (4.17) may suggest. 

The cl protein is found in numbers typically less than 180-350 for a healthy cell, and 
error from counting statistics alone contributes zby^N) to the intrinsic noise. The uncer- 
tainty in AG and the dimerization equation of state also induces a spread of trajectories. A 
convenient way to catch these and other sources of stochasticity is to assert a diffusion ten- 
sor, Dij, with i,i twice indexing over the protein species. This is equivalent to considering 
some Langevin-like equation but allows us to use partial differential equations to explore 
ensembles of trajectories. Noise of greater complexity then Gaussian noise can be readily 
studied with the software developed. 

As with the other parameters in this model they are left determined by experiments 
and have been given a simple form to understand the effects that chance events could 
have on the dynamical picture. In particular, it can wash away higher resolution features 
which can, for example, remove the gap between two nearby stabilities, facilitating perhaps 
the appearance of bistability when the dynamical picture may suggest more structure. 
From the modeling perspective it can be useful to use large, perhaps unphysically large 
(analagous to a high temperature limit), diffusion to speed up simulations and see some of 
the more basic underlying features of the phase space. As diffusion gets smaller, the steady 
states take exponentially longer to converge and can quickly become computationally 
prohibitive. In these scenarios, where converging to a steady state from any arbitrary initial 
condition is computationally prohibitive, the steady states from larger diffusion cases can 
serve as useful initial conditions since they have already undergone a large degree of 
this convergence process and the remaining time to equilibrium is for the finer resolution 
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features to develop since the coarse features were already present at the beginning of the 
simulation. 

4.4 Little's Mutants 

Mutations in the Or region have been performed, see Little et al. [1999], and we would 
like to check the performance of our model and suggest means to calibrate accordingly. 
The qualitative features compare to experiment but quantitative comparison is premature; 
in particular, the time and concentration units need to be understood better. As is, the first 
passage times predict a smaller barrier then observed. 



Diffusion Tensor 


mutant 


Dxx 


^yy 


Xi^iWild-type 


3.90778 


2.22668 


Ai2l 


0.775855 


2.22668 


A323 


2.6462 


1.5024 


^^323' 


2.0016 


2.86862 



Table 4.1: Values for diffusion tensors, in Little et al. [1999] is nearly A222/ the cl 

affinities are identical and Cro values differ by -.6 kcal/ mol in the first and third site. 



4.4.1 Wild-Type Evolution 

Four trajectories of distributions are depicted in Tables A.5, A.6, A.7 and A.8, corresponding 
to the wild-type and mutants A121, A323, and A323/. It is evident that in the first 200 seconds, 
an initially peaked distribution has developed a second peak and after about 600 seconds 
the distributions have reached their steady states. 

4.4.2 Steady state properties for Or Mutations 

The Tables A.9, A.IO, A.ll, A.12, A.13, A.14, A.15, A.16, A.17, A.18, A.19, A.20, and A.21 
archive the results. Depicted are the phase space portraits alongside the Probability and 
Flux Landscapes and above the passage time distributions to each of the three stationary 
points. Assorted expectation values are calculated and can be compared with the wild-type. 

In particular, we notice that with diffusion the peaks of the solutions are considerably 
displaced towards smaller numbers from where one would expect the peaks from the 
deterministic equations. This is required to satisfy the boundary condition, and it would be 
interesting to see how a Langevin-based approach can reproduce this ensemble property 
in that different trajectories would require correlation. 

Noise in this model seems to destabilize at least the mutants Am, A131, and A312; either 
no gap, or a small gap is present. The phage seems most robust in response to mutations 
in Ori in that A*23 all have a sizeable gap. 
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4.5 Hybrid particle and ensemble approaches 

It is important to see how different paths contribute differently to the overall energy and 
entropy profiles, see Seifert [2005], Stock et al. [2008]. To illustrate the shape of these 
landscapes it is sufficient to evaluate the steady state ensembles values along determinstic 
trajectories. These calculations are depicted in Tables A.25, A.26, and A.27. ^ It is evident 
that different paths can converge to different locations and demonstrates the motivation 
for a (more fundamental) dynamic partition function. 

4.6 Variations on a theme 

We revisit the equations of motion with RNAp still included in the partition function. Here, 
the promotional states are the natural ones, the binding of RNAp to prm for cl promotion 
and the binding of RNAp to pr for Cro promotion. These promotional states enter the 
dynamics as follows. 

djcl) _ 1 . 
dt 

mpi * Trm * Eci{[RNAp][cI]e-^ + [RNAp] [c7] + [RNAp][cI][Cro]e-^) 

ry ^'-^33 ^^103 ^^203 ^'^023 

+ Trmu^Eci{[RNApYe~^ + [RNAp][cI]e ^ + [RNAp][Cro]{e ^ +e ^] 



AG223 



[RNAp][cI][Cro]e-^ + [RNAp][CroYe-^ + [RNApje'^ 



(4.15) 



djCro) _1 
dt ~7S 

AGgg ^'-^31 ^'^32 o 

mpi * Trr * Ecro * {\KNAp\e~^ + [RNAp] \cl\e~^ + [RNAp] \Cro\e~^ + {KNApY 
^ Cro 

^Cro 

(4.16) 

The Z in these equations is now from 2.7, with values listed in table 4.2. The inclusion 
of these changes leads to much higher numbers for cl, up to 1500 proteins, though the 
lysogen peak still dominates. The time scale is slowed down in comparison to the dynamics 
without RNAp. The \KNAp\ = 600nM was chosen to be in the range of the experiment 

^Seifert [2005] suggests a definition of path entropy as s(t) = — Zn((x(T), r)), with x(t) some 
stochastic trajectory. The definition is sufficient to define a path entropy for any function x(t), 
regardless of the motivation of such functions. The principle of maximum caliber, see e.g. Stock 
et al. [2008], is an effort to define a dynamic partition function, over all such paths. There has been 
some success with this method for low number of states, but in the discrete case, or the continuous 
case here, the implications are largely unknown. 
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Klumpp and Hwa [2008] and is presumed to be a fixed number for the simulations. The 
results are displayed in tables A.22 and A.23. 

An important contribution of the codes written in this thesis is the capacity for noise 
with spatial and temporal noise. As an example, a spatially dependant noise that distorts 
locally the classical kinetics is given in 4.17. A multiplier of 1000 is used to control the 
numerics. ^ 



^Computationally, diffusion equations with finite elements requires a very fine grid, or a rotating 
coordinate system when the noise is smaller or comparable to the underlying flow. The factor of 
a thousand is an artificial enhancement to control the numerics, (i.e. a high temperature limit of 
sorts) In general, it is best to fit all parameters simultaneously with experiments. 




(4.17) 



iz£ 
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Ai23/ '^ild — type Oj^ affinities (kcal/mol) 


AG,-yfc 


protein 


Ori 


Or2 


Or3 


ref. 


equation 


monomial 


premotes 


AGooi/AGoio/AGioo 


-15 


-13 


-12 


Zhu et al. [2004] 


AGqqq = 0.0 


1 




AG002/AG020/AG200 


-18.4 


-17.1 


-19.5 


Zhu et al. [2004] 


AGqqi = -12.0 


in] 




AGcoop — —6.9 


AGqiq = -13.0 






pnc = 2.5 * 10~^ 


AGiQQ = -15.0 


[r/l 




[RNAp] = 300,600nM Klumpp and Hwa [2008] 


AGqq2 = -19.5 


\Cro] 




AG30 = —11.5, —10.5, RNAp binding to Cro promoter. 


AGq2Q = -17.1 


\Cro] 




AG03 = —11.5, —13.5, RNAp binding to cl promoter. 


AG2QQ = -18.4 


\Cro] 
















r ^ -1 9 

[CtY 






1 A n 

^^01 — ^^100 + ^^001 








AG022 — AG020 + AG002 


[CroY 






/\ ( -riT^r* /\ I -Tr»'^r» —\— /\ I -Tnr»r» 

^^220 — ^^^020 ^ ^^^200 


[Lroj 






AG202 — AG20O ~^ AG002 


[Lroj 






AG120 — AGioo + AG020 








AG210 — AG200 ~^ AGqio 








AG102 — AGioo + AG002 


fry] rrrni 






AG201 — AG200 ~^ AGqoi 


fry] ICrn] 






AG012 = AGqio + AG002 


\Ci] \Cro\ 






AG021 — AG020 ~^ AGqoi 


\Ci] \Cro] 

L J L J 






A A 1 A 1 A 1 A 
ALrlll = ALrlQQ + ALrQlQ + ALrQQl + l^^cOOp 


\Ci\'^ 






A A 1 A 1 A 

ALr222 = ^^200 + ^^020 + ^^002 








A A 1 A 1 A 1 A 

AG112 = AGiQQ + AGqiq + AGcoo^ + ^^002 








AG121 = AGiQQ + AGq2Q + AGqqi 


[LZJ [Lroj 






A A 1 A 1 A 1 A 

ALr211 = A^200 + A^OOl + A^OIO + ^^coop 


[L,zj [uroj 






AG22I = AG2QQ + AGq2Q + AGqqi 


[Lzj [Lroj 






AG212 = AG2QQ + AGqiq + AGqq2 


[CzJ [Lro\ 






AG122 = AGiQQ + AGq2Q + AGqq2 


[CzJ[CroJ^ 






AG3Q = -10.5 


1 






AGq3 = -13.5 


1 






AG33 = AG3Q + AGq3 


1 






AG31 = AG3Q + AGqqi 


[cl] 






AG32 = AG3Q + AGqq2 


[Cro] 


Zv 




AG1Q3 = AGlQQ + AGq3 


[cl] 


Zx 




AG2Q3 = AG2QQ + AGq3 


[Cro] 






AGqi3 = AGqiq + AGq3 


[cl] 


Zx 




AGq23 = AGq2Q + AGq3 


[Cro] 


Zx 




AG113 = AGiQQ + AGqiq + AGq3 


[clY 


Zx 




AG123 = AGiQQ + AGq2q + AGq3 


[cl] [Cro] 


Zx 




AG213 = AG2QQ + AGqiq + AGq3 


[cl][Cro] 


Zx 




AG223 = AG2QQ + AGq2Q + AGq3 


[CroY 


Zx 



Table 4.2: Values of variables referenced by the partition function with RNA polymerase 
states included. A.22 uses the values of [RNAp],AG3o, and AG03 listed first, whereas A.23 
and A.24 use the values listed second. 



5 Conclusions and Future Directions 



Dafi es mir - oder Allen - so scheint, darausfolgt nicht, dafi es so ist. 
Ludwig Wittgenstein On Certainty 

We have compared the probability that an unknown scientist should have 
found out what has been vainly sought for so long, with the probability that 
there is one madman the more on the Earth, and the latter has appeared to us 
the greater. 

Henri Poincare Science and Hypothesis 

The particular contributions of this author towards this thesis, and towards 
other dynamical systems not discussed in this thesis, are archived in Borggren 
[2011]. This has allowed for much deeper quantitative inspection of the dynamical 
systems with noise then hitherto discussed with the phage A. This has allowed 
for the prediction of behaviour of phenotypes not yet synthesized and further 
comparison with those that have been. Some subtlety involved with the removal 
of RNA polymerase from the partition function has been discussed. The A phage 
has been discussed in physical and computational terms in hopes to present the 
system to a physics, mathematics, and computational audience. The A phage is 
proving to be to biophysics what hydrogen was for quantum mechanics ^ and as 
such is a perfect testing ground for empirical and ab initio physical, computational 
and mathematical methods for gene expression. 

Efforts should be made to link the connection with the decision portion of the 
cell cycle to the bigger picture involving the upstream and downstream reactions 
and understand the connections between all the parameters and variables, see 
St. Pierre [2009], Arkin and Ross [1998]. Solutions are only initial conditions for 
the next portion of the cell cycle and our simple models are in danger of being 
too simple, as we are not left with much room under the rug to sweep degrees of 
freedom we choose to ignore. 

In particular, volume should be connected into the modeling with greater care. 
The binding and unbinding are in the end very specific bond types occurring at 
very specific times, and proceeding in terms of idealized concentrations tends to 
wash out these subtleties. I cannot say with confidence if the model here can or 

^ A point Jin Wang has made to motivate this study. 
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cannot corroborate the findings of St. Pierre [2009] with respect to the sensitivity 
of the lysis and lysogeny decision being governed by volume. It can be said 
though, that the decision being sensitive to volume is not mutually exclusive to 
the decision also being controlled by protein concentrations. In the mathematical 
sense, we expect the spaces to be connected; volume is in the end the denominator 
of these densities. Much of the fluctuation can be considered fluctuation of specific 
volume of cl. For example, if one were to carefully wrap a cylinder around the 
helix-bend-helix motif of the protein and associated dimerizations and tabulate 
the volume, that volume would vary greatly over the course of, say, a binding and 
unbinding process. This model lacks any specific protein configuration and in 
general, specific volume can vary a great deal over this space. Discrete jumps in 
that space would appear as noise in this model and would need to be tabulated as 
such in the diffusion tensors, for example. From a phenomenological perspective 
it is a useful construct that is parameterized, ideally, by direct measurement and is 
a highly intuitive way to reproduce the general features of the circuit. 

5.0.1 On Computation and Information 

It is important to realize that computers evolve simultaneously in two directions. 
There is the physical interface, the its that make the bits and the its that control 
them, but there is also the human interface that establishes how we control the its 
that control the bits, and how the consequent bits are presented back to us in such 
a way usable to our central nervous system. 

5.0.1.1 The Physical Interface 

The physical interface evolves in the direction of greater complexity. Moore's law 
illustrates famously that rates of (serial and electronic) computation are limited by 
size and heat production. Parallelization, reversible logic, and quantum computa- 
tion have all been under continual development to circumvent these issues and 
with considerable success. Where a decade ago one would speak of how fast their 
computer was, now one speaks of how many cores it has. This is parallelization's 
contribution. Hq core would be the limiting case if this trend of the processor to 
duo — core to quad — core to 2^ core persists, but appended to the technical diffi- 
culties involved with space are now the technical difficulties involved with time, 
the synchronization of the bits from all the its. The limitation that won't go away 
here is the speed of light; if A needs a signal from B before time t, the situation 
is hopeless unless B is spatially located within ct of A. Multiple processors on a 
chip enable a reduction in heat dissipated by providing channels for which heat 
can be carried away from a processor. Cells to cells exhibit analogous, although 
much more sophisticated signaling, are far closer to the Hq — core limit, and being 
waterproof would certainly be convenient for transporting heat elsewhere; 98.6 
degrees is certainly quite a clue that bodies are quite adept at managing this 
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Figure 5.1: Varia- 
tional forms used in 
this work. 



Python code 



a = v*u*dx + 0.5*k*(inner(grad(v) , -velocity*u)*dx + 

inner (grad(v) , D*grad(u) )*dx) 
L = v*u0*dx - 0.5*k*(inner(grad(v) , -velocity*u0)*dx + 

inner (grad(v) , D*grad(u0) )*dx) 



housekeeping heat. Reversible logic evades the heat issue while exacerbating the 
spatial issue since intermediate steps are kept. 

5.0.1.2 The Human Interface 

The human interface evolves in the direction of greatest simplicity. Quantum Com- 
puter, classical computer, topological quantum computer, chemical computer,... is 
trivia far from the concern of the user who ultimately just wants to push a but- 
ton and get a reliable answer. Contrast the Fortran punch codes of old with the 
following lines of python code that defines our variational forms: 

These lines of code are the same no matter what mesh we choose, velocity we 
choose, time step, initial condition, and diffusion tensor. When the code is called 
we better have defined what they mean, but until then they are freely transf err able. 
See more examples in Logg et al. [2012]. 

A Gui was built to provide buttons to push to generate dynamical trajectories, 
Langevin trajectories, mesh trajectories (in which an entire domain is propagated 
in accordance with the flow field), and diffusion trajectories. Systems built are the 
A phage, A phage with RNA polymerase, A phage approximating Oi states, the 
cyclin system (see Wang et al. [2008]), the Lorenz attractor, and can operate with 
imaginary numbers for the ode case. 

Here is a screenshot showing multiple instances of the Gui from Borggren 
[2011]. 

Figure 5.2: This 
is the ^^00000 Gui 
demonstrating 
user interaction 
with equations 
of motion. It 
reads equations of 
motion of a system 
and builds a Gui 
to interact with 
that system 
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5.1 On the importance of generality and the abstract 

It is important that methods and calculations be reproducible and transferrable to 
different systems. The study of dynamics in the most general settings cannot be 
neglected for proper understanding of the circuitry and signaling of cells. 

Sun Ra says that it is laughable when knowledge is attributed to a man, 
Ra and Mugge [1980], and this rings all the more true when the discussion 
is one that seeks to synthesize elements from the great discoveries of physics: 
quantum mechanics and statistical mechanics, with the great discoveries of biology, 
evolution, genes, and DNA to just begin. Chemistry has found itself central 
to the discussions of biology as a result of being best positioned and adept at 
comprehending the molecular basis and even the purest of mathematics at this 
stage cannot escape its insights from being applied to the problems of life. At 
this interface; engineers, medical doctors, biologists, physicists, mathematicians, 
chemists, logicians, theologians and everybody else with hands and a mind are 
trying to contribute to this story. The vacuousness of the borders of our disciplines 
has become evident and we find Nature ready and willing to adapt to any of our 
biases. The wisdom of Poincare on mathematics equally applies to the whole of 
science as well: 

We know very well that mathematics will continue to develop, but 
we have to find out in what direction. I shall be told "in all directions," 
and that is partly true; but if it were altogether true, it would become 
somewhat alarming. Our riches would soon become embarrassing, and 
their accumulation would soon produce a mass just as impenetrable as 
the unknown truth was to the ignorant. 

So it is the case that Science and Method begins with a selection of facts, see 
Poincare [2010]. For Poincare this consisted in those with the greatest chance of 
recurrence and to each of us our lives come to be mirrors of our particular criterion 
of this selection process. Wittgenstein in his youth would further require that this 
selection be from all the facts, see Wittgenstein [1922], but by the end of his term 
he would have brushed that off to some of the symbolic games that we play, see 
Wittgenstein [1975]. The world is more akin to a collection of statements then a 
collection of facts and a statement can be true, false, or meaningless. Perhaps one 
may want to continue this line of thinking to get a world that is a collection of 
words, but words require letters and something to denote nothing to tell where 
one word ends and the next word begins. The statistical mechanics of random 
words and strings of letters may be equally important, see Stephens and Bialek 
[2008, 2010]. 

What sets the subjects apart on this stage is no longer the content of their 
subject matter - it is one and the same - but the psychological predispositions of 
the practitioners. What are the relevant degrees of freedom? No two agree, but 
nor is it necessary for them to do so. A mutual respect is in order, and as the 
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tricks of noise and dynamics that can crash markets are now being applied to the 
problems of life, let us walk gently into this night, none of these days have been 
promised to us. Recall a giant meteor from outerspace was required to wipe out 
the dinosaurs and open up our niche. Nobody can look a three-headed pterodactyl 
with flourescent wings in the eyes and call themself the pinnacle of creation, and as 
far as I can tell, their is nothing prohibiting DNA to make such a genome. Science 
is still new, our ignorance will be infinite, even while our knowledge can be infinite 
too. For want of a common goal let us choose, at least, a moral one. 

The only thing that could unite the planet is a united space program 
[health center]... the earth becomes a space station [health center] and war 
is simply out, irrelevant, flatly insane in a context of research centers, space- 
ports, and the exhilaration of working with people you like and respect toward 
an agreed-upon objective, an objective from which all workers will gain. Hap- 
piness is a by-product of function. The planetary space station [health center] 
will give all participants an opportunity to function. 

William S. Burroughs, The Place of Dead Roads [NAB] 

Toss in some clean water, money, and the cure for pain and we may asymptoti- 
cally approach a civil society. When the subject matter is infection and disease and 
the consequences are pain and joy it would be imprudent for physics to wait for a 
topological quantum computer and all orders of precision worked out to chime 
in to this story, which is still caught up in tabulating everything. From the looks 
of the White House, and the dazzling noise of stock markets apparently stuck 
between radio stations, the stalemate seem to persist straight up to the top. 

The challenges of this century can provide an incomparable opportunity to- 
wards global scientific collaboration as we learn to harness the resources we are 
enveloped in, or we can watch ourselves fumble the oil into sea, and another 
mysterious earthquake, tornado and tsunami wipe out our cities, as we choose the 
right color for our data points, and argue over what data points to delete. 

The human genome is at least the magnitude of mystery it was a decade 
ago and bioinformatics has been more successful at deciphering the code than 
electricity, magnetism, statistical mechanics, quantum mechanics, and all the other 
tools physicists pridefully declare universality with. Some basic assumptions 
clearly need to be relinquished. 

5.2 On the importance of the particular and the measurable; Case Studies 

A model of the phage A has been analyzed to illustrate how a particular system 
is queried. Given that these crystal structures have been measured and good 
empirical force fields exist, these values of AGs can in principle begin to be 
predicted and compared to the values deduced from experiments. Ab initio 
calculations of the Or region from first principle Hamiltonians may take some time 
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but until then it is exactly these types of mechanisms that motivate experiments 
like NSLS to do such crystallographic measurements that empirical protein and 
DNA force fields are developed to illuminate. Understanding the relation of 
these nanoscopic structures and dynamics to the macroscopic flows of energy and 
entropy is the basic question necessary to address in order to appreciate the full 
value of the accumulated data and begin to use it for clinically relevant purposes. 
Nuclear magnetic resonance, NMR, is also regularly used to investigate more 
directly these ensemble properties in vivo, near equilibrium, and can also be used 
to constrain these numbers. The difference between in vitro and in vivo values for 
the weights has already been demonstrated to induce topological changes. These 
numbers have a taxonomy tree of their own as measurement techniques get refined. 
Our phenomenology here continues to rely on fits to the primary datum from the 
biological measurements. Physics will learn a lot from the networks that make up 
biological systems and the manners in which they give rise to emergent degrees of 
freedom. In order to address the concerns of Little et al. [1999], and perhaps many 
other scientists who are typically more concerned about large scale qualitative 
differences, like the difference of lysogeny and lysis, we will scan parameters 
as much as we can in this volume of spacetime. The particular values of AG 
are dependent in some complicated way to evolution, volume, temperature, pH, 
doping, in vitro, in vivo and a host of other unknowns, can and do change on 
large scales. SchrodingeraAZs anthropomorphic remark can be applied here: it 
is exactly robustness to these types of changes that allow it to exist at 300K for 
so long in the first place. Here we note a difference in the use of the terminology 
aAYrobustaAZ between biologists and the use in subsequent theoretical work. 
For Little, robustness is the property that two states, lytic and lysogenic, are still 
present after large scale mutations, changes in AG due to switching, say, site 3 
with a copy of site 1. This is a transformation on the level of the partition function 
in that the AG values change. This differs from the use of robustness in the sense 
of Zhu et al. [2004], Wang et al. [2010], Aurell and Sneppen [2002], which uses this 
word to describe the small rate at which a cell state switches spontaneously as a 
result of noise. The phage can be robust in both senses. It is even robust against 
our attempts to model it. 

Binding studies for drug docking often involve phenomena of this sort and 
computational drug design as well. This is a largely trial and error process, and 
perhaps some topological information is due. Biologists are often skeptical of the 
latest grand unification scheme from the physicists and the applied mathematicians, 
and their skepticism is mimed here. The main point of separation, as I see it, 
between the two subjects is in a sense a question of topology. The mass of an 
electron is a number with small error bars that with hard work can be made even 
smaller. It is repeatable and exists on other planets too, the numbers of biology are 
scarce when present. The facts they prefer to select on are topological in nature: to 
lyse or not to lyse? 
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5.2.1 Robustness, robustness, stability, and stability 

I am referring to something that mere words will never be able to 
express, relative, absolute, full, empty, still alive and no longer alive, because, 
sir, in case you dont know it, words move, they change from one day to the 
next, they are as unstable as shadows, are themselves shadows,... 
Jose Saramago, Death with interruptions 

In the overlap of subjects it is quite expected for there to be ambiguity between 
the different uses of the same word. This is all too commonplace and has been a 
perennial source of confusion in this study. Two examples are worthy of note, the 
notion robust in the title of Little et al. [1999] is different than the use of the same 
word in the subsequent work of Wang et al. [2008], Zhu et al. [2004]. The first use is 
to denote the propensity of phage A to maintain a working switch after particular 
mutations in the sites. This is a change at the level of the partition function; the 
genome has fundamentally changed, yet the bistability persists, it is robust. The 
second use refers to something fundamentally different: the capacity of a lysogen 
to stay a lysogen by continuous maintenance by cl protein to keep, in particular, 
the cro gene from expression. The use of the word stability in Little et al. [1999] is 
exactly what Wang et al. [2008], Zhu et al. [2004] refer to as robustness, reserving 
their use of the word stability to refer to the mathematical use of describing the 
nature of the stationary points. Stability and robustness are used interchangeably 
in Aurell and Sneppen [2002]. 

With a four order magnitude spread in the refereed values of the numbers 
quantifying robustness, it is possible that there is some misunderstanding here. 
Indeed the numbers reported as experimental in Wang et al. [2010] are not strictly 
speaking experimental numbers at all, they involve simulation as well, and Little 
has now decreased those numbers to 10~^ or 10~^, in units of lysogens flipping per 
generation, but has not published these numbers in a refereed journal and caveats 
the measurement as very difficult when it arises. Finding a more realistic diffusion 
tensor and better oligomerization equations of state are likely to help better 
quantify the robustness in the model. We note that in the large concentration limit, 
where the dynamics picture is sufficient, the phase space completely partitions 
into two basins of attraction and the number is simply 0. The stability puzzle is 
still puzzling. 

5.3 Algebraic, Analytical, or Computational Methods? 

A computational approach, the finite elements, were eventually decided upon for 
the practical solution of the differential equations herein Logg et al. [2012]. Other 
methods were available and attempted. In particular, an algebraic approach was 
also investigated (see Marcus [1960], Figueiredo et al. [1998]), and a brute-force 
analytic approach was used in Snider [1999]. These are briefly discussed but were 
found not as immediately illuminating as the numerical approach throughout. 
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5.3.1 N on- Associative Commutative Algebra 

A differential equation is associated with a non-associative algebra such that in- 
variant manifolds of the differential equation correspond to invariant sub-algebras 
of the algebra. In particular, given a dynamical system, expand in a form, 

^i = L^ijUh''' (5-1) 

j=l k=l 

we can then define a matrix M = BA and a system Uj = UfMijUj. This is then a 
quadratic system and the methods of Marcus apply; we can define an algebra by 
introducing a product Uj * Uj = ^i^ik^kj + ^ij^jk)- 

5.3.2 Recursion Equations for Expansion coefficients 

In particular, a dynamical system with some noise is concisely written adopting 
the summation convention 

dtP = diDijdjP-di{FiP), (5.2) 

Roman indices run from 1 to K over the generalized coordinates of the system. This 
relates the time change of the surface P with a diffusion term and the flow field of 
the dynamical system. The noise has effectively been smoothed and encoded in 
the surface P through the diffusion tensor D/y. To be explicit, by 3/D/y9yP we mean 
MDijidjP)). 

The zero is used to denote the time index and thus K + 1 indices over the 
non-negative integers are summed over. Ansatz P, a polynomial expansion of the 
form 

p(xo,^i,..vxk) = pzo,ii,...,/K n ^^-^^ 

and similarly for the components of F 

fm(^0/^l/-/^K) =/m,/o,/i,...,/x Yl ^) ' (^-4) 

We have allowed the index corresponding to powers of x/ to run over a countably 
infinite set if need be, however the m counts only the K indices of the coordinates. 
In general we need an analogue expansion for Dmn- We will soon specify to the 
constant diffusion case, which is lowest order of the general case 

^mni^O/ ^1/ "'/ ^k) — ^m,n,io,ii,...,iK 1 1 "^/* i^-^) 

Term by term we will decipher equation (5.2), 
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dtP{xo,xi,...,XK) = Pio,i^ i^{iox'° Yl xl = pi^,+ij^ /k(^o + 1) n 

i=\,...,K i=0,...,K 

(5.6) 

Here we specify constant diffusion tensors. For the next term in equation (5.2), the 
coefficients on the diagonal are related by 



;=i ^ 

and off diagonal by 



+i)(f„+i) n (5-8) 

The last term in equation (5.2) has two contributions through the product rule. 



(dniFm)P{xQ,Xi,...,XK) = ifmj'^j'^ 4+1 /^(4 + l) Yl ^j)Vk,h im k H ^) ' 

i=0,...,K j=l,.-,K 

(5.9) 



■1— r ^ -1— r i- 

FmdmP{xo.Xi,...,XK) = ifm/^/^ i^^ Yl ^pPkA + Yl ' 

(5.10) 

These become 



{dnrFm)P{XQ,Xi,...,XK) = fnr,i'i'/^^^^ Y\ (^-H) 

;=i K 



FrndmP{xo,Xi,...,XK) = f„rj'i'^,.,j'^Pi^^^^^^^^ Yl xj '. (5.12) 

Let us simplify by saying the diffusion tensor is diagonal. Substituting into the 
differential equation yields 
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Pio+l,h /k(^0 + 1) n ^? = ^rnmPk,h i„r+2 i^iim + 2){im + 1) fl ^'j 

j=0,...,K i=l,:;K 

i'+i" 

f' ...A+i,...,i^pf"f"..x,..v,i;^(4 + 1) n ^j^l^f) 

;=i K 

- fm£j'i'Pi''i'{ i,n+V',...,i'liC + '^) n ^/ '■ 

M K 

Setting the coefficients of Ylj=o,...,K equal on both sides gives the recursion 
relation 



Pio+Ui,-,iK (^'O + 1) = DmmPio,^,...,iM-,iK + 2) (im + 1) 

~ E /m,/[,/i,...A+l,...,/5.P/[,'//,...,4, /^(4 + l) (5-14) 

ij+ij -1, 



A steady state solution in particular requires = for all i. So that the index can 
be dropped altogether. 



= DmmPh,...,i„+2,..JK{im + 2){im + 1) 

- E fm,i[ 4+1 /jcP/;',...,4,...„/5^(4 + 1) (5-15) 



I] fm,i[ i'i.Pi'{,...,im+V',..J^ (4 + 1) • 



The software of Borggren [2011], can also perform calculations based on these 
approaches. 

5 A From biology to physics to biology 

Physics has always relied heavily on progress from biology for its own progress. 
The lung predates the steam engine and it is unlikely the latter would have 
emerged without the scientific inspection of the former. The central nervous 
system predates the understanding of electricity and again the latter would not 
have emerged without the scrutiny of the former. 

Consider the timelime of the study of metabolism, see Kilgour [1955]. In around 
1500, Leonardo da Vinci likened the burning of a candle to animal nutrition, from 
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which he inferred that an atmosphere that could not support combustion would 
also not support animals. In 1648, Van Helmont coined the word 'gas' in a 
journal for medicine under the title Ortus Medicinae. The first gas described 
as such was CO2, the biproduct of animal respiration. From 1660-1678 Robert 
Boyle developed experiments that studied the changes in a volume of air due to 
respiration and combustion and similarly for Mayow's Treatise on Respiration. In 
1754 Black discovered fixed air, CO2. In 1766, Cavendish discovered inflammable air, 
H2. In 1772, Rutherford described residual air, N2. In 1780, with their new invention, 
the calorimeter, Lavoisier and Laplace concluded that da Vinci's hypothesis of 
respiration being a form of combustion is indeed the case in their beautiful paper 
Memoir on Heat. Thus the groundwork had been laid for thermodynamics and the 
subsequent industrial revolution, the driving force of which having been the study 
of respiration. 

One only need to read the bibliography of Hermann Helmholtz to see that 
physiology was a driving force in the study of electricity and acoustics, see 
Helmholtz [1954]. 

Progress in physics in the 20th century, given the equipment and good math- 
ematical abstractions of the 19th century, managed on its own for a century 
without much aid from biology and chemistry, but now it is time to revisit the old 
phenomena, with our new lights, and actually solve some problems. 

Biologists and chemists have handed us the structures of their molecules, the 
wirings of their systems, and the phase diagrams of their compounds. If those 
are the answers and statistical mechanics and quantum mechanics provide the 
questions, what is left but the proof? 
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Appendix; Tabulation of Results 



Ineluctable modality of the visible: at least that if no more, thought through 
my eyes. Signatures of all things I am here to read, seaspawn and seawrack, 
the nearing tide, that rusty boot. Snotgreen, bluesilver, rust: coloured signs. 
Limits of the diaphane. But he adds: in bodies. Then he was aware of them 
bodies before of them coloured. How? By knocking his sconce against them, 
sure. Go easy. Bald he was and a millionaire, maestro di color che sanno. 
Limit of the diaphane in. Why in? Diaphane, adiaphane. If you can put your 
five fingers through it, it is a gate, if not a door. Shut your eyes and see. 
James Joyce, Ulysses 

We live, I regret to say, in an age of surfaces. 
Oscar Wilde, The importance of being earnest 

A.l Results 

The phase portraits for the 27 mutations are shown in Tables A.l, A.2 and A.3. 

The switch component of 3 of the 9 mutations in Table A.l are destabilized: 
Aii2,Ai22/ and A132. Transients however are still largely important for the out of 
equilibrium switch properties. 3 of the 9 mutations in Table A.2 are destabilized: 
^212/^221/ ai^d A231. Transients however are still largely important for the out of 
equilibrium switch properties. 5 of the 9 mutations in Table A.3 are destabilized: 
^311/^321/ ^322/ A331 and A332. Transients however are still largely important for the 
out of equilibrium switch properties. 

The locations of the stationary points are tabulated in Tables A.l, A.2, A.3, and 
A.4. 

Time evolution of distributions is given in Tables A.5, A.6, A.7 and A.8, corre- 
sponding to the wild-type and mutants A121, A323, and A323/. 

The Tables A.9, A.IO, A.ll, A.12, A.13, A.14, A.15, A.16, A.17, A.18, A.19, A.20, 
and A.21 archive the results of mutants with three stationary points. Depicted are 
the phase space portraits alongside the Probability and Flux Landscapes and above 
the passage time distributions to each of the three stationary points. Assorted 
expectation values are calculated and can be compared with the wild-type. 

Representative trajectories, encoding global features, are depicted in Tables 
A.25, A.26, and A.27. 
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Figure A.l: Con- 
tours with station- 
ary Ci and Cro pro- 
duction with Or-^ 
fixed. The axis 
are in protein num- 
bers. 3 of the 9 mu- 
tations are desta- 
bilized: Aii2,Ai22a 
and Ai32 
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Figure A.2: Con- 
tours with station- 
ary Ci and Cro 
production with 
Or^ fixed. The 
axis are in protein 
numbers. 3 of the 
9 mutations are 
destabilized:A2i2,A22i, 
and A231. 
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Figure A.3: Con- 
tours with sta- 
tionary Ci and 
Cro production 
with Or^ mutated 
to Or^. 5 of 
the 9 mutations 
are destabilized: 

A311,A32l, A322, A331 

and A332. 
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Properties of the Stationary Point 


Mutant 


X (Ci proteins) 


y (Cro proteins) 


A112 




404.90 


0.0009 


A122 




404.89 


0.0232 


Al32 




404.89 


0.1173 


A212 




329.88 


0.02954 


A22I 




0.3535 


712.24 


'^231 




1.7522 


312.15 


A311 




3.4415 


214.70 


A32I 




1.7532 


312.14 


-^322 




10.253 


309.57 


'^331 




6.1379 


153.09 


A332 




20.403 


152.35 



Table A.l: Mutants with one stationary point. 



Properties of Stationary Points 


Mutant 


index 


X (Ci Proteins) 


y (Cro Proteins) 


A113 







644.05 


0.00023 




1 


0.31299 


310.32500 



Table A.2: Mutants with two stationary points. 
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Properties of Stationary Points 


Mutant index x (Ci Proteins) y (Cro Proteins) 


Am 







5.6833 


301.287 




1 


1.8180 


310.218 




2 


127.87 


0.03220 


A121 







19.9875 


332.348 




1 


0.82427 


462.571 




2 


127.86 


0.8235 


A123 wild-type 







0.141753 


462.57273 




1 


644.05 


0.00611 




2 


29.979 


178.229 


-^131 







30.0039 


191.3784 




1 


3.43968 


216.27658 




2 


127.812 


4.168042 


^^133 







48.6464 


130.6941 




1 


0.6369 


216.311 




2 


644.043 


0.0309 


A211 







17.505 


206.876 




1 


37.3015 


4.8364 




2 


0.8244 


462.570 


^213 







0.14175 


462.572 




1 


595.05 


0.00679 




2 


29.938 


177.711 


A222 







37.234 


590.517 




1 


2.6864 


712.240 




2 


329.884 


0.7555 



Table A.3: Mutants with three stationary points. 
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Stationary Points 


Mutant 


index x (Ci Proteins) 


y (Cro Proteins) 


A223 







0.06002 


712.247 




1 


595.04 


0.17382 




2 


83.624 


103.62 


A232 







56.762 


289.227 




1 


10.272 


312.072 




2 


329.829 


3.82207 


^^233 







0.30934 


312.153 




1 


146.0568 


74.57208 




2 


594.8508 


0.87984 


A312 







19.08756 


102.4368 




1 


13.3744 


163.4361 




2 


158.8459 


0.47835 


A313 







44.4330 


100.41068 




1 


0.63744 


216.3115 




2 


436.73205 


0.054279 


^^323 







0.30935 


312.15312 




1 


133.92 


52.23903 




2 


436.115 


1.393207 


A333 







1.24472 


153.1449 




1 


267.603 


28.44186 




2 


417.693 


7.90558 



Table A.4: Mutants with three stationary points (continued) 
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time 

t = l 


+0 

2.76s^0^^^^M94^^^^^nS9 0.264 0.379 


rime Evolution 
+100 

J.a2^^^^a0W2^^^^50KM 0.015B 0.0211 


+200 

260^^^^0W3^^^^0.0065i 0.00987 0.0132 


f = 301 


V 

}.63^^^^0m37^^^W0.00750 0.0113 0.0150 


V 

2.6e^^^^^^3a^^^Wq.00770 0.0116 0.0154 


V 

2.25^^^^aoo3a^^^2i.oo7d7 0.01 15 o.orsj 


t = 601 


V 

7.69^^^^a^^7^^^^^m75& 0.0114 0.0152 


V 

2.2s^^^^aoro^^^^^oo74a 0.0112 0.0150 


V 

4.9S^^^^^^36^^^^^073e 0.0111 0.0148 


f = 901 


V 

a.23^^^^^^^^^^^^^0730 0.0146 


V 

MI^^^^^^^^^^^^^07M 0.0145 


V 

f.52^^^^^^^^^^^^g007M ftOfOe o.or44 



Table A.5: The time evolution of a gaussian initial condition relaxing into the steady 
state for the wild type. The higher peak corresponds to lytic while the wider peak is 
lysogenic. 
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time 

t = l 


+0 

■l.52e-137 ft Ml 0.2S1 0.422 0.562 


rime Evolution 
+100 

■1.236-105 0.0120 0.0241 0.0361 0.0481 


+200 

8.506-120 0.00824 0.0165 0.0247 0.0330 


i = 301 


}.36s-lU 0.00783 0.0157 0.W35 0.0313 


8.63S-109 0.00974 0.0195 0.0292 0.0390 


3.946-108 0.011} 0.0223 0.0334 0.0446 


t = 601 


266s-m 0.0119 0.0239 0.0358 0.0477 


9.886-109 0.0123 0.0247 0.0370 0.0493 


3.436-109 0.0125 0.0250 0.0376 0.050} 


i = 901 


].41s^^^^^^^^^^^^0.0252 0.0378 0.0503 


7.56e-U0 0.0126 0.0252 0.0378 0.0503 


5.246^^^^^^^^^^^^0.025} 0.0377 0.0502 



Table A.6: The time evolution of a gaussian initial condition relaxing into the steady 
state for the 121 mutant. The dominant peak is the lysogenic peak. 
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Time Evolution 


time 

t = l 


2.5Be-213 0.116 


0.236 0.355 0.473 


+ 1UU 

4.3BS-96 0.00775 0.0155 0.0232 0.0310 


+ZUU 

Jl 

l.lls-69 0.00600 0.0120 0.0160 0.0240 


t = 301 


J 

2.54S-61 0.00835 


L 

0.0167 0.0251 0.0334 


4.52S-5S 0.00944 0.01S9 0.02S3 0.037S 


9.74S-57 0.00998 0.0200 0.0299 0.0399 


t = 601 


J 

3.55S-S6 0.0103 


0.0205 0.030e 0.0410 


6.59S-56 0.0104 0.0209 0.0313 0.0417 


Jl 

216S-55 0.0105 0.0211 0.0316 0.0422 


f = 901 


J 

4.5e^^^^^O010^^^ 


0.02 r3 0.0319 00425 


7.66^^^^^^^^^^^^ 0.0214 0.0321 0.0428 


Jl 

r.os^^^^^^^^^^^^ ao2r5 0.0323 0.0431 



Table A. 7: The time evolution of a gaussian initial condition relaxing into the steady 
state for the 323 mutant. The dominant peak is lytic. The gap is reduce compared to 
wild type. 
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time 

t = l 


+0 

9.746-209 0.106 0.217 0.325 0.433 


Time Evolution 
+100 

7.95S-83 0.00655 0.0131 0.0197 0.0262 


+200 

2.63S-60 0.00422 0.00845 0.0127 0.0169 


i = 301 


JL 

9.24S-54 0.00401 0.00601 0.0120 0.0160 


2.21S-5 1 0.00462 0.00923 0.0138 0.0185 


1.75S-50 0.00496 0.00993 0.0149 0.0199 


t = 601 


JL 

3.64S-50 0.00516 0.0103 0.0155 0.0207 


JL 

5.70S-50 0.00529 0.0106 0.0159 0.0211 


JL 

7.92S-50 0.00535 0.0107 0.0161 0.0214 


i = 901 


JL 

i.oe^^^^^^^^^^^^o.oioe 0.0162 0.0216 


JL 

1.34^^^^^^^^^^^^0.0108 0.0163 0.0217 


JL 

1.60^^^^^^^^^^^^^.0109 0.0163 0.0217 



Table A.8: The time evolution of a gaussian initial condition relaxing into the steady 
state for the 323p mutant. The dominant peak is lysogenic. 
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A121 steady state and first passage time distributions 



Ori 



TATCACCGCCAGAGGTA 



Or2 



TAACACCGTGCGTGTTG 



Or3 



TATCACCGCCAGAGGTA 



{1,2,1} 




-AG 



Ori Or2 Or3 



100 150 200 250 300 350 





A121 Assorted Expectation Values 



J^xp{x,y)dxdf ~ 
J^p{x,y)dxdy 



ratio to Ai; 



< X > (Ci protein) 



''^yp{x,y)dxdy 
]^p{x,y)dxdy 



114.8' 



1.234 



< y > (Cro protein) 



!ia,Mx,y)pi.^,y)dxdy 

J^p{x,y)dxdy 
J^Ti{x,y)p{x,y)dxdy 

J^p{x,y)dxdy 
J^T2{x,y)p{x,y)dxdy 

Jf^ p{x,y)dxdy 



193.112 



1.137 



< To > (seconds to ts) 
< Ti > (seconds to lytic) 
< T2 > (seconds to lysogen) 



11177.264 
103986.923 
22020.823 



0.527 
1.113 
0.083 



Table A.9: Properties of mutant A121 
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Ai23 Steady state and first passage time distributions 



Ori 



TATCACCGCCAGAGGTA 



TAACACCGTGCGTGTTG 



Or3 



TATCACCGCAAGGGATA 



^{1,2,31 




Ori Or2 Or3 






5 Assorted Expectation Values 



Jp xp{x,y)dxdy/ J^p{x,y)dxdy 



< X > (Ci protein) 



j'^TQ{x,y)p{x,y)dxdy/ j'^p{x,y)dxdy 
Sp Ti{x,y)p{x,y)dxdy / fpp{x,y)dxdy 
In y)p(.^' y)dxdy/ p{x, y)dxdy 



93.076 



< y > (Cro protein) 



< To > (seconds to lytic) 



93387.22 



< Ti > (seconds to lysogenic) 
< T2 > (seconds to transition state) 



263858.749 
21214.802 



Table A.IO: Properties of wild-type A123 
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Ai3i steady state and first passage time distributions 



Ori 



TATCACCGCCAGAGGTA 



Or2 



TATCACCGCAAGGGATA 



Or3 

^{1,3,1} 



TATCACCGCCAGAGGTA 



-AG 




Ori Or2 Or3 






Ai3i Assorted Expectation Values 



ratio to Xy. 



J^p{x,y)dxdy 
]^p{x,y)dxdy 

J^p{x,y)dxdy 
J^Ti(x,y)p{x,y)dxdy 

j^p{x,y)dxdy 
J^T2{x,y)p{x,y)dxdy 

J^p{x,y)dxdy 



< X > (Ci protein) 



117.454 



1.262 



< y > (Cro protein) 



143.183 



0.843 



< To > (seconds to ts) 
< Ti > (seconds to lysogenic) 
< T2 > (seconds to lytic) 



8084.495 
13499.860 
12908.285 



0.381 
0.051 
0.138 



Table A.ll: Properties of mutant A131 
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Ai33 Assorted Expectation Values 


ratio to A123 


< X > (Ci protein) 


J^xp{x,y)dxdy 
J^p{x,y)dxdy 


93.227 


1.002 


< y > (Cro protein) 


j^p{x,y)dxdy 


137.624 


0.8106 


< To > (seconds to ts) 


in'^o{x,y)p{x,y)dxdy 
J^p{x,y)dxdy 


12264.638 


0.578 


< Ti > (seconds to lytic) 


J^Ti{x,y)p{x,y)dxdy 
J^p{x,y)dxdy 


21628.286 


0.232 


< T2 > (seconds to lysogenic) 


J^T2{x,y)p{x,y)dxdy 
J^p{x,y)dxdy 


198518.768 


0.752 



Table A. 12: Properties of mutant A133 
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A211 steady state and first passage time distributions 



Ori 



TAACACCGTGCGTGTTG 



Or2 



TATCACCGCCAGAGGTA 



Or3 



TATCACCGCCAGAGGTA 




Ori Or2 Or3 



100 150 200 250 300 





A211 Assorted Expectation Values 



J^xp{x,y)dxdf ~ 
J^p{x,y)dxdy 



ratio to Ai; 



< X > (Ci protein) 



''^yp{x,y)dxdy 
]^p{x,y)dxdy 



329.760 



3.543 



< y > (Cro protein) 



!ia,Mx,y)pi.^,y)dxdy 

J^p{x,y)dxdy 
J^Ti{x,y)p{x,y)dxdy 

J^p{x,y)dxdy 
J^T2{x,y)p{x,y)dxdy 

Jf^ p{x,y)dxdy 



120.570 



0.710 



< To > (seconds to lytic) 
< Ti > (seconds to ts) 
< T2 > (seconds to lysogenic) 



7388.126 
39165.111 
122448.847 



0.079 
1.846 
0.464 



Table A.13: Properties of mutant A211 
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A213 steady state and first passage time distributions 



Ori 



TAACACCGTGCGTGTTG 



Or2 



TATCACCGCCAGAGGTA 



Or3 



TATCACCGCAAGGGATA 






3 Assorted Expectation Values 



J^xp{x,y)dxdf ~ 
J^p{x,y)dxdy 



ratio to Ai; 



< X > (Ci protein) 



''^yp{x,y)dxdy 
]^p{x,y)dxdy 



105.183 



1.130 



< y > (Cro protein) 



!ia,Mx,y)pi.^,y)dxdy 

J^p{x,y)dxdy 
J^Ti{x,y)p{x,y)dxdy 

J^p{x,y)dxdy 
J^T2{x,y)p{x,y)dxdy 

Jf^ p{x,y)dxdy 



141.037 



0.831 



< To > (seconds to lytic) 
< Ti > (seconds to lysogenic) 
< T2 > (seconds to ts) 



90915.013 
212601.812 
20741.420 



0.974 
0.806 
0.978 



Table A. 14: Properties of mutant A213 
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Table A. 15: Properties of mutant A222 
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Table A.16: Properties of mutant A223 
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Table A. 17: Properties of mutant A232 
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Table A.18: Properties of mutant A233 
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Table A.19: Properties of mutant A312 
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Table A.20: Properties of mutant A313 
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Table A.21: Properties of mutant A323 
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Time Evolution 


time 


+0 


+1200 


+2400 


f = 


i 

X 


-5.a3^^^^^^^^^^^^?.57e-05 2.36S-05 3. !4e-05 


A 

-3.(3^^^^^^^^^^^ J.Me-OS 177S-0S 2.36S-0S 








t = 3600 


-4.77^^^^^^^^^^^^^^^ l.SOs-05 2.00S-05 


-l.7a^^^^^^^^^^^^99e-0« r.35e-05 ?.aOe-05 


-3.4I^^^^^^^^^^^J.36e-0« h25s-0S I.67S-05 








t = 7200 


■4.59^^^^^^S^^^^^U^-05 1.73S-0S 2.306-05 


-5.04^^^^^^^^^^^^5e-O5 2.329-05 3. lOs-05 


-4.S'0^^^^^R37^^^^^r.«7e-O5 2.S1S-0S 3.75S-05 








t = 10800 


u 


u 


u 


■4.43^^^^^^^^^^^^^3s-05 IPe-05 4.25S-05 


-3.M^^^^^^^^^^^^3re-05 3.47s -05 4.62s -05 


-3.23^^^^^^^^^^^ 2.44S-05 3.67S-05 4.S9S-05 



Table A.22: The time evolution of a gaussian initial condition relaxing into the steady 
state for the wild type with RNA polymerase in the partition function and constant 
diffusion tensor. Notice that the two states, lytic and lysogenic, are gapped only for 
the window of time from 4800 to 8400 seconds. Although the emphasis in the literature 
is with steady states, it is an interesting possibility that the two state phenomena is a 
strictly out of equilibrium phenomena, existing only for a window of time in the cell. 
A bacteria cell cycle is around 1800 seconds, suggesting that for many distributions the 
steady state is not reachable in the necessary time scale. 
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Time Evolution 


time 


1 n 


+ lzUU 


+Z4UU 


t = 1200 








-3.5^^^^^^^^^^^^^.24e-(M 7.S6S-06 L05e-05 


-l.0Si^^^^^^^^^^^^4.67e-(M 7.01S-06 9.34S-06 


















t = 4800 








■L09^^^^^^^^^^^^^e-O6 6.25e-06 e.34e-06 


■7.62^^^^^^^^^^^^V5s-06 6.US-06 SABs-06 


-4.S'7^^^^^^^^^^^^.aite-(M 5.95s -06 B.lOs-06 
















t = 8400 








I. I5^^^^^^^^^^^^.46e -C« 5. 76s -06 07s -06 


■1.70^^^^^^^^^^^^. 16s -06 5.59S-06 S.Ols-06 


■2.04^^^^^^^^^^^^^^ -06 6.2 J s -06 S.97S-06 
















t = 12000 


■2.23S-06 a.l4s-07 3.B6e-06 6.90S-06 9.95S-06 


■2.33^^^^^^^^^^^^^Ss-06 7.2SS-06 l05s-05 


■2.3SS-06 9.06S-07 4.I9e-06 7.4SS-06 J.Oas-05 









Table A.23: The time evolution of a gaussian initial condition relaxing into the steady 
state for the wild type with RNA polymerase in the partition function and constant dif- 
fusion tensor with different parameters from A.22 . The lysogenic peak in the classical 
dynamics is near 1500 cl proteins. To find a distribution with a second peak required 
this enlarged domain and was facilitated by an initial conditition with many cl (3500). 
The short axis is 0-1000 Cro proteins, and the long axis is 0-5000 cl proteins. One can 
see the global effects through the large displacement of the lysogenic peak from its 
stationary point. 
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time 

f = 


+0 

1 


Time Evolution 
+600 

-3.49S-27 3.329-06 9.979-06 I33e-0S 


+1200 

-5.29^^^^^^9^6^^^^n9-06 }. 029-05 }. 369-05 


i = 301 


- }. S5^^^^^U3^^^^S. S6s -06 1.33S-05 }.77s-05 


- 1.02^^^^^^e^^^^l329-05 1.979-05 2.639-05 


-4. 1 r^^^^^^^^^^^^gte -05 3.309-05 4.409 -05 


t = 601 


-2.49^^^^^^^^^^^ 4.71S-05 7.07S-05 9.43S-05 


u 

-4.36^^^^^^^^^^^^000}U 0.000 W5 0.000233 


-3.26^^^^^^^^^^^^^fi^^^^^m403 



Table A.24: The time evolution of a gaussian initial condition relaxing into the steady 
state for the wild type with concentration-space dependant diffusion. The peak corre- 
sponds to the lytic state. Two strong tails develop, one extending towards potentially 
another peak. 
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Table A.25: Entropy and Energy densities evaluated along example paths for mutations 
in Or2,r3 with the Ori fixed. Different trajectories can have different limiting values of 
entropy and energy density depending on the path taken. 
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Table A.26: Entropy and Energy densities evaluated along example paths for mutations 
in Or2^r3 with the Ori mutated into Or2- Different trajectories can have different limit- 
ing values of entropy and energy density depending on the path taken. 
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A3*, Path Entropy and Path Energy Densities 
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Table A.27: Entropy and Energy densities evaluated along example paths for mutations 
in Or2^r3 with the Ori mutated into Ors. Different trajectories can have different limit- 
ing values of entropy and energy density depending on the path taken. 



Sometimes the glass unbreaks and falls to heaven. 
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